reflecting the partitioning of shared food resources. Dietary protein exhibited the theoretically predicted negative allometric relationship with body mass, whereas digestible-energy content scaled positively. Our results demonstrate cryptic size-based separation along spatial and dietary niche axes-despite superficial similarities among species-consistent with the idea that body-size differentiation is driven by selection for divergent resource-acquisition strategies, which in turn underpin coexistence. Foraging and space-use behaviors were nonlinearly related to body mass, supporting the hypothesis that behavior scales allometrically with size. However, explaining the variable functional forms of these relationships is a challenge for future research.

Behavioral variability is a major determinant of ecological performance and fitness, both within and among species (Relyea, 2001). For example, behavioral plasticity enables animals to adjust their resource use in response to competitive pressure (Miner et al., 2005). However, the range of possible behaviors is constrained by functional traits, especially body size, which sets physical limits on what animals can do (Bonner, 2011;Peters, 1983;Schmidt-Nielsen, 1984). Competitive pressures should therefore select for body-size differentiation among sympatric species and play an eco-evolutionary role in shaping species coexistence and community assembly (Anaya-Rojas et al., 2021;Dayan & Simberloff, 2005;Pelletier et al., 2009). Still, there are many groups of sympatric species that differ in size yet superficially appear to use more-or-less the same resources. One explanation for this phenomenon is that size-dependent morphological and physiological traits lead to subtle differences in behavior that collectively differentiate species' niches (du Toit & Olff, 2014;Jarman, 1974). Inconspicuous relationships between body size and behavior may thus be central to understanding the organization of animal communities.

Many morphological and physiological traits scale allometrically (nonlinearly) with size (Capellini et al., 2010;Peters, 1983;Schmidt-Nielsen, 1984). These relationships are described by the exponential function y = a  x b , where y is the trait, x is body mass, a is a proportionality coefficient, and the scaling exponent b ≠ 1 (as opposed to isometric, or linear, scaling where b = 1). Because resource acquisition and risk avoidance are linked to traits that scale allometrically, such as metabolic rate and locomotion, it has been proposed that behaviors should exhibit similar scaling relationships (Dial et al., 2008). If so, then allometries of behavior may also underpin niche differentiation and coexistence. Yet, evidence for allometric scaling of behavior remains scarce (Cloyed et al., 2021;Laca et al., 2010;Preisser & Orrock, 2012;Sensenig et al., 2010).

Several factors make it difficult to evaluate the functional form and biological significance of associations between body size and behavior. Unlike many morphophysiological relationships, behavioral scaling is expected to be inherently noisy due to factors other than body size that affect animals' decisions and performance (Dial et al., 2008). One such source of noise is local environmental context (Cloyed et al., 2021). Another is phylogenetic constraint (Capellini et al., 2010;Shine, 1994). In ungulates, for example, differences among lineages in life history, gut anatomy, dentition, social structure, and other traits complicate the investigation of scaling relationships (Clauss et al., 2013;Demment & Van Soest, 1985;Illius & Gordon, 1992). These potentially confounding influences can be controlled by comparing closely related sympatric species (Grant & Grant, 2008;Losos, 2009;Shine, 1989), thereby minimizing environmental and phylogenetic variability and honing the focus on body size per se as a driver of behavior. An added benefit of this approach is that it can explicitly link allometries of behavior with niche differences among co-occurring populations, thereby shedding light on mechanisms that promote body-size differentiation and stabilize coexistence in assemblages of closely related species.

Yet another challenge in evaluating behavioral allometries is to separate the influences of opposing ecological forces that might amplify the noise in (or at least complicate the interpretation of) any observed relationships. Two such factors are resource acquisition and predator avoidance, which interact to determine space use in most animals (Brown & Kotler, 2004). In ungulates, as in many other taxa, vulnerability to predation is inversely related to body size (Sinclair et al., 2003). Moreover, predator-avoidance strategies vary qualitatively with body size: many small species rely on crypsis and hiding, many mid-sized species rely on gregariousness and flight, and adult megaherbivores are essentially invulnerable to predators (Atkins et al., 2019;Dial et al., 2008;Ford et al., 2014;Jarman, 1974;le Roux et al., 2018). Such variability in risk-avoidance behavior might obscure allometric scaling in behaviors driven by nutritional requirements, or vice versa, and to the extent that these two forces drive behavior in similar directions, it would be difficult to attribute causality to either one. Thus, it is useful to study behavioral allometry in ways that can isolate the effects of fear and hunger. One way to do this is to focus on systems with weak predation pressure-which are increasingly common worldwide as carnivore populations decline (Ripple et al., 2014)-where the necessity of risk avoidance is reduced and habitat selection should be more unambiguously resource driven. Experimental manipulations of habitat attributes associated with risk avoidance and/or resource availability can be used to test this assumption.

African savanna ungulates are an ideal system for studying behavioral allometry and its influence on niche differentiation, for multiple reasons. First, assemblages are diverse and contain closely related sympatric species that differ in size, which has long drawn the attention of niche theorists (Hutchinson, 1959). Second, the behaviors of African ungulates are well characterized, at least in broad strokes (Estes, 2012). Third, allometric scaling is hypothesized to influence foraging behavior and dietary niches in these species. For example, the Jarman-Bell principle (Bell, 1971;Jarman, 1974) posits an inverse relationship between body size and diet quality; this relationship is often observed (Clauss et al., 2013;Kleynhans et al., 2011;Owen-Smith, 1988;Potter et al., 2022;Sensenig et al., 2010) and is thought to arise from allometries of food intake, mass-specific metabolic requirements, and/or other traits associated with feeding and digestion (Demment & Van Soest, 1985;Müller et al., 2013;Potter & Pringle, 2022;Steuer et al., 2014). Fourth, savannas exhibit both coarse-and fine-scale spatial heterogeneity in the quantity and nutritional quality of the plants that ungulates eat (Cromsigt & Olff, 2006;Hopcraft et al., 2010), which stems from interactions among geology, rainfall, and fire (Pringle et al., 2016;Scholes, 1990;Smit et al., 2013) along with the ecosystem-engineering effects of mound-building termites (Davies, Baldeck, et al., 2016;Okullo & Moe, 2012). Allometric scaling of behavior could thus emerge from interactions among size-varying traits that directly affect animal morphophysiology and habitat characteristics, such as resource distribution, that contextualize the relationship between animals and their environment (Dial et al., 2008;Milne et al., 1992;Peters & Wassenberg, 1983).

We investigated how body size influences space use and foraging by three sympatric, congeneric species of spiral-horned antelopes (Tragelaphus spp.)-bushbuck (T. sylvaticus), nyala (T. angasii), and greater kudu (T. strepsiceros)-in Gorongosa National Park, Mozambique (Figure 1). These closely related browsers (most recent common ancestor 6-7 million years ago; Hassanin et al., 2018) have essentially nonoverlapping body-size distributions, spanning a four-fold range in average adult female body mass (bushbuck roughly half the size of nyala, nyala roughly half the size of kudu) and differing up to six-fold among individuals (Kingdon, 2015). All three species co-occur in southeastern Africa, occupy similar habitats, eat similar plants, and are thought to compete (Coates & Downs, 2005;Ehlers Smith et al., 2020;Fay & Greeff, 1999;Tello & Van Gelder, 1975). Gorongosa historically supported robust populations of bushbuck, nyala, and kudu (Tinley, 1977), but they and other large mammals were nearly extirpated during the Mozambican Civil War (1977-1992). Since that time, populations have increased; as of 2018, the three focal species occurred at roughly equal densities (~1 km À2 ) in our study area at the core of the park (Stalmans et al., 2019).

Two features of Gorongosa, both shared by many African savannas, are key to our study design. First, apex predators were scarce during our study, and Tragelaphus spp. experienced little predation (Methods: Study system), enabling us to focus on the allometry of resource-acquisition behaviors while minimizing the potentially confounding effects of risk. Second, fungus-farming termites (Macrotermitinae) create marked spatial heterogeneity in resource availability (Figure 2; Appendix S1: Figure S1). Termites concentrate moisture and nutrients in the soils around their nests (termitaria, "mounds"), leading to distinctive plant assemblages that are productive and nutrient rich relative to the surrounding matrix (Joseph et al., 2013;Seymour et al., 2014;Sileshi et al., 2010). Mature mounds are typically separated by tens to hundreds of meters (Pringle & Tarnita, 2017) and may thus serve as "islands" of palatable food for ungulates (Grant & Scholes, 2006;Holdo & McDowell, 2004;Levick et al., 2010;Okullo et al., 2013). However, mounds are not always heavily used by herbivores (Davies, Levick, et al., 2016;Muvengwi et al., 2013Muvengwi et al., , 2019;;Van der Plas et al., 2013), and factors other than food-such as risk avoidance or microclimate-might provide alternative explanations for high herbivore activity around mounds (Anderson et al., 2016;Joseph et al., 2016).

We collected hourly GPS locations from adult females of each species in each of 2 years. We paired these spatial data with (i) habitat classifications derived from remotely sensed imagery, (ii) surveys of woody vegetation, (iii) diet-composition data from fecal DNA metabarcoding, (iv) nutritional analyses of plants and diets, and (v) experimental manipulations of forage availability to test three general hypotheses and nine specific predictions (Table 1).

These hypotheses are motivated and unified by the framework outlined above; collectively, they predict that allometries of behavior shape patterns of resource use and partitioning in heterogeneous landscapes, with the sign of allometric relationships depending on the behavior. For example, we hypothesized that all antelope species are attracted to termite mounds by the localized availability of high-quality forage, but that the strength of selection for mounds scales negatively with body size, because small animals can (and may need to) subsist on high-quality diets, whereas larger individuals can (and may need to) tolerate lower quality diets (Bell, 1971;Clauss et al., 2013;Jarman, 1974;Potter & Pringle, 2022). Conversely, home-range sizes and step lengths should scale positively with body size, because localized resource hotspots are sufficient to fulfill the dietary requirements of small animals, whereas larger individuals require more food and must range farther to get it (Harestad & Bunnell, 1979;Illius & Gordon, 1987;Noonan et al., 2020). For similar reasons, we expected the scaling of several behaviors to differ between wet and dry seasons: larger antelopes should relax selection for mounds in the wet season, when forage is abundant in the matrix. Finally, we hypothesized that size-dependent differences in space use are associated with differences in diet composition and the partitioning of shared food resources. We predicted that all three species predominantly eat mound-associated plants, but that selectivity for these taxa is weaker in larger animals, which should eat more plant species but have less energy-and protein-rich diets. Thus, bushbuck and kudu should have the most dissimilar diets and nyala should be intermediate. Support for these predictions would indicate that the scaling of behavior with body size promotes separation along spatial and dietary axes in ways that stabilize the coexistence of closely related species (Table 1).

We conducted fieldwork from 2014 to 2016 in the south-central portion of Gorongosa National Park, Mozambique (Figure 1a). Mean annual precipitation is 850 mm, most of which falls in the wet season from November to March (Tinley, 1977). Annual rainfall during our study was 1200 mm in 2014, 688 mm in 2015, and 754 mm in 2016 (mean 881 mm). The study area, defined by the movements of GPS-collared antelopes (Figure 1b), is bounded on the northeast side by Lake Urema and includes part of the Urema floodplain. From north to south, the floodplain grades into seasonally flooded savanna dominated by fever trees (Acacia syn. Vachellia xanthophloea) and palms (Hyphaene coriacea), and then into woodlands (Acacia-Combretum savanna, sand forest) with patches of saline grassland (Daskin et al., 2016). Termite mounds created by Macrotermes mossambicus and M. subhyalinus are a conspicuous and abundant feature of the landscape, covering 1.1% of the study area (6.9 of 629 km 2 ) with a mean density of 68 km À2 (Figure 1c,d). These roughly conical termitaria (which can exceed 5-m height and 20-m diameter) support dense woody thickets (Tinley, 1977), including trees up to 25-m tall (Figure 2; Appendix S1: Figure S1). Mature termite mounds are spatially overdispersed at local scales (mean nearest-neighbor distance ~50 m; Appendix S1: Figure S2) but aggregated at very large scales (Tarnita et al., 2017).

Gorongosa's large-mammal populations declined by >90% during the Mozambican Civil War (Stalmans et al., 2019). By 2016, herbivore biomass had recovered to nearly prewar levels, but with altered relative abundances. Buffalo (Syncerus caffer), hippopotamus (Hippopotamus amphibius), and elephant (Loxodonta africana) dominated the prewar assemblage, whereas mid-sized ungulates accounted for the majority of biomass during our study (Stalmans et al., 2019;Tinley, 1977). Predator recovery has been slower. Of four top carnivores present in the 1970s, F I G U R E 2 Spiral-horned antelopes (Tragelaphus spp.) on a termite mound in Gorongosa. This composite image, created from a series of camera-trap photographs, shows females of all three focal species (bushbuck at bottom and top; nyala at bottom left, center, and right; kudu at center left) browsing on the woody vegetation characteristic of Macrotermes spp. mounds (image courtesy of Jennifer A. Guyton). Additional images illustrating termite-induced heterogeneity in Gorongosa are given in Appendix S1: Figure S1.

T A B L E 1 Hypotheses and predictions tested in this study.

Specific predictions Support Evidence S3; Tables S6 andS7) only lion (Panthera leo) persisted through the war. Some 65-80 lions were known to be alive during our study (Bouley et al., 2018;Pringle, 2017), less than half the prewar population size (Tinley, 1977), and these ranged across our study area (Figure 1e). Leopard (P. pardus) and African wild dog (Lycaon pictus) were extirpated during the war and were not present during our study, but have subsequently been reintroduced, starting in 2018 (Bouley et al., 2021).

Our focal antelope species are woodland-affiliated, ruminant browsers that rely on concealment for predator avoidance (Estes, 2012). Bushbuck are solitary or paired; nyala and kudu occur in herds of 3-12. All three species overlap in space and time in our study area (Figure 1b). The adult female weights measured in this study (mean AE SD: bushbuck 33 AE 4 kg, nyala 73 AE 9 kg, kudu 120 AE 14 kg; Figure 1f) are similar to those reported elsewhere (Kingdon, 2015). In the absence of leopard and wild dog, and with human hunting curtailed (Bouley et al., 2018), lion were the only major potential predator of Tragelaphus spp. during our study. However, lion diets from 2012 to 2020 (n = 307 kills; Bouley et al., 2018Bouley et al., , 2021) ) were dominated by warthog and waterbuck (79% of kills), did not include bushbuck (no kills), and only rarely included nyala (seven kills, 2%) or kudu (three kills, 1%). We thus infer that the real per capita risk of predation for these species was negligible during our study (although this does not necessarily preclude perceived risk and associated avoidance behaviors, which may be "hard-wired" to some extent; Berger et al., 2001).

To test the hypothesis that termitaria support distinctive woody-plant assemblages with greater abundance, richness, and nutritional quality than the matrix, we conducted vegetation surveys in June-July 2015 and June 2016. These data also enabled us to quantify the relative availability of different food-plant taxa, which we used to analyze selectivity. We selected 30 points along roads within the minimum convex polygon (MCP) bounding all antelope GPS locations (Figure 1b). Points were spaced evenly along each road (550-1700 m apart) and spanned the full spectrum of vegetation densities (0%-92% canopy cover) and fire-return intervals (1.3-17 years from 2000 to 2015, as per Daskin et al., 2016) in the MCP. From each point, we walked perpendicularly to the road on both sides for a randomly selected distance between 10 and 250 m and found the nearest termite mound. From the center of that mound, we walked in a randomly selected direction and distance between 20 and 50 m to locate the center of a paired matrix plot.

Mound plots (n = 60) encompassed the entire surface area of the mound. We walked around the mound edge, where the difference in topography was apparent (Appendix S1: Figures S1 andS2), and recorded the circumference (from which we calculated radius) with a GPS. Mound height was estimated (to nearest m) for all but two; for those, we used the mean height of the other 58. We estimated mound area as the lateral surface area of a cone:

where h and r are mound height and radius, respectively (area range 25-539 m 2 , mean AE SD = 177 AE 109 m 2 ; cumulative area surveyed = 10,750 m 2 ). Matrix plots (n = 60) were circular with 8-m radius (plot area = 201 m 2 ; cumulative area surveyed = 12,060 m 2 ) to approximate the mean surface area of mound plots.

In each plot, we censused the overstory (woody plants ≥2-m tall). We identified plants using keys (Coates Palgrave, 2002;van Wyk, 2013) in consultation with botanists. Uncertain identifications were recorded using "cf." to denote similarity with a known taxon or "morpho." for morphospecific labels. We measured basal area of all stems at 20-cm height. We photographed the canopy at the center of each plot using a 15-mm lens on a cropped-sensor digital SLR camera (24-mm full-frame equivalent) pointed straight up at a height of 1 m. We then quantified canopy cover in ImageJ software (Schneider et al., 2012) by converting images to grayscale and using a brightness threshold to classify pixels as plant or sky. We used Welch's unequal-variance t-tests to analyze the effect of habitat type (mound vs. matrix) on canopy cover, species density (species per unit area), and basal-area density across the 120 plots (the latter two variables were square-root transformed for normality). To more fully account for differences in plot area and stem density between mounds and matrix, we also compared species richness using individual-based rarefaction (at 342 individuals, the number sampled in matrix plots) in EstimateS v9.1 (Colwell & Elsensohn, 2014).

We investigated variation in foliar-nutrient concentrations for nine common woody-plant species. We collected leaves from individuals growing on and off mounds, in areas of higher and lower fire frequencies, along each of three roads in the MCP. We tried to collect paired mound and matrix samples of each species in each fire frequency along each road, but this was not always possible. In total, we sampled 76 individuals (range 4-12 per species), with each species represented from both mounds (range 1-6) and matrix (range 3-6) and high (range 2-6) and low (range 1-6) fire frequencies (except for Combretum imberbe, which we did not find on mounds). We sampled green leaves at heights accessible to all antelope species (≤1.5 m). Leaves were dried at 50 C and analyzed at the Cornell University Nutrient Analysis Laboratory (Ithaca, NY) for %N, C:N, B, Na, Mg, Al, P, S, K, Ca, Mn, Fe, Ni, Cu, Zn, and Ba. We tested whether nutrient concentrations differed depending on (i) whether a plant species was mound versus matrix affiliated (based on where the majority of records occurred in the vegetation surveys; Appendix S1: Table S1) and local growth conditions, including (ii) habitat type (mound vs. matrix), (iii) fire frequency (high vs. low), and (iv) road identity (a proxy for spatial heterogeneity in factors that might affect soil and plant nutrient contents, such as flood regime). After a significant multiple analysis of variance (MANOVA) on all 16 nutrients (square-root transformed) with these four factors (Wilks' Λ = 0.14, F 80,269 = 1.70, p = 0.0009), we analyzed each nutrient in a separate analysis of variance (ANOVA) with the same factors (n = 76 measurements each).

In June 2014 and July-August 2015, we chemically immobilized adult females of each antelope species by darting them using species-specific combinations of thiafentanil, medetomidine, and azaperone. All procedures were approved by Princeton University's Institutional Animal Care and Use Committee (protocol 1958-13) and conformed to guidelines from the American Society of Mammalogists (Sikes, 2016). We weighed each bushbuck and nyala (nearest 0.1 kg) and estimated kudu weight as a function of chest girth, based on data from the other two species. We collected a fecal sample and fit each animal with an Advanced Telemetry Systems G2110E iridium GPS collar. In all, we collected data from 57 individuals: 19 bushbuck (11 in 2014, 8 in 2015), 16 nyala (10 in 2014, 6 in 2015), and 22 kudu (12 in 2014, 10 in 2015). Collars recorded locations hourly and transmitted data daily to a server via satellite. Mean GPS measurement error for these collars at our site is ~13 m (Atkins et al., 2019). Collars were remotely released when they entered low-battery status (usually 10-12 months after deployment) and retrieved without recapturing animals. To minimize error in habitat-selection analyses, we followed Lewis et al. (2007) and Long et al. (2014) in excluding GPS locations that had both a two-dimensional fix and a dilution of precision >5 (which affected <1% of all fixes). Another 0.14% of fixes were removed as erroneous because they fell well outside the park boundary.

We assessed habitat selection relative to three factors: distance to nearest termite mound, woody cover, and lion utilization. We used a continuous metric of mound proximity in lieu of a categorical on/off mound variable to reduce bias from GPS error and also because the effects of termite mounds on soils and plants typically extend well into the matrix (Baker et al., 2020;Pringle et al., 2010;Sileshi et al., 2010) such that herbivores may use mound-associated resources even when not on a mound (although the spatial extent of this effect in our system is unknown). We included the latter two covariates to help control for residual variation in selection for termite mounds, given the known effects of tree cover on real/perceived predation risk and antelope behavior (Atkins et al., 2019;Ford et al., 2014;le Roux et al., 2018;Tambling et al., 2013;Valeix et al., 2009). Including lion utilization also helped to test our assumption that antelope movements were driven primarily by resource distribution rather than predation risk from the sole extant large carnivore.

We mapped mound distribution using airborne Light Detection and Ranging (LiDAR; Davies et al., 2014;Levick et al., 2010) data collected in August 2019 by Wooding Geospatial Solutions (Everton, South Africa). Flights were conducted within 2 h of solar noon at 880 m above ground level, yielding terrain-elevation measurements at 50-cm resolution. Using LiDAR-derived digital terrain models and the hillshade tool in ArcGIS, we manually digitized termitaria locations and sizes based on differences in slope and shape (Appendix S1: Figure S2). Although we did not quantify the accuracy of this approach, a previous study that used automated classification to map termitaria in similar habitat from LiDAR data with coarser resolution detected 78%-90% of mounds >0.5-m tall (Davies et al., 2014), which are those likely to be used by antelopes.

To map woody cover, we used a supervised classification of 1.8-m resolution satellite imagery (WorldView-2, Digital Globe, Longmont, CO) collected in July-August 2010 to categorize each pixel as either woody (overstory) or herbaceous (understory) vegetation. The resulting layer was accurate in comparison with a visual classification of 300 randomly selected points (accuracy 87%, sensitivity to woody cover 79%, specificity 92%; Appendix S1: Table S2).

To map lion utilization, we used locations of GPS-collared lions in Gorongosa (Bouley et al., 2018) to estimate 100% fixed-kernel utilization distributions (UDs) with a 216-m resolution (smoothing factor set to 60% of the reference bandwidth; Kernohan et al., 2001). UDs quantify the relative intensity of space use by an animal within its home range (Millspaugh et al., 2006). We used data from eight collared females, including one from each of the five prides that consistently occupied the MCP in 2015 (the middle of our study) and three from a sixth pride in which animals ranged widely and often independently of each other both inside and outside the MCP. Following Valeix et al. (2009) and Davies et al. (2016a), we averaged the 100% UDs across individuals to produce a population-level UD.

To quantify antelope habitat selection, we estimated resource selection functions by fitting generalized linear mixed models (GLMMs; Bolker et al., 2009;Gillies et al., 2006;Zuur et al., 2009) with a binomial error distribution and logit link function to used (antelope GPS points, coded "1") and random (available habitat, coded "0") locations in a use-availability design (Johnson et al., 2006;Manly et al., 2002). We used the near function in ArcGIS 10.0 to calculate the distance between each used location and the nearest termite mound. Next, we spatially joined antelope GPS locations to the lion UD layer in ArcGIS and computed proportional woody cover in a fixed radius (20 m for bushbuck and nyala and 270 m for kudu, based on the analysis in Appendix S1: Table S3) around each antelope location using the "extract" function in the raster package (Hijmans & van Etten, 2014) in R v3.3.1 (R Core Team, 2015). Random locations were spatially joined to each habitat layer as described above for antelope locations.

Because habitat selection is a scale-dependent process (Johnson, 1980), we quantified selection for mounds by each antelope at each of two spatial scales: (i) the area encompassed by the 629-km 2 MCP around all antelope GPS locations ("landscape scale" below), and (ii) the home range of each individual ("home-range scale" below, also known as third-order selection). We estimated 95% fixed-kernel home ranges using the adehabitatHR package in R to quantify home-range-scale selection. We followed Long et al. (2014) in determining the number of random locations required to adequately represent available habitat at each scale.

We fit the binomial GLMMs for each species-scale combination in lme4 (Bates et al., 2015). We standardized and centered all predictor variables by subtracting their mean and dividing by their standard deviation (Cade, 2015;Kutner et al., 2004), which placed them on the same scale and enabled direct comparison of effect sizes. We verified that no pair of predictors exhibited problematic collinearity (jrj ≤ 0.40 for all variables). Our approach was similar to that of Long et al. (2014). We began by fitting a global model that included fixed effects for all predictors (Mound, distance to mound; Woody, proportional woody cover; Lion, lion utilization), along with a random intercept grouped by individual ID (to account for serial autocorrelation in each animal's GPS location data). For simplicity, this approach assumes a linear relationship between each covariate and the probability of use by antelope, although we acknowledge the possibility of nonlinear functional forms. We then determined the optimal random-effects structure for each combination of species and scale by comparing (using AIC c ) the random-intercept model to three additional models, each of which included a random slope for one of the three predictor variables (Zuur et al., 2009). In all cases, inclusion of random slopes greatly reduced AIC c relative to the random-intercept model (ΔAIC c > 1936; Appendix S1: Table S4), indicating substantial variation in habitat selection among individuals. Therefore, we included random slopes for each predictor in subsequent models of habitat selection for each species and scale.

To account for temporal variation in habitat selection (Burkepile et al., 2013;Spitz et al., 2018), we investigated interactions between each predictor and two categorical variables: Season (dry, April-October; wet, November-March) and Time (day, 5:00 AM to 5:00 PM; night, 5:00 PM to 5:00 AM). We then used AIC c to compare fully parameterized models from the first stage of analysis (i.e., those including all three predictor variables as fixed effects and random slopes grouped by individual) with two additional models for each species and scale: one including all pairwise interactions between the three original predictors and Season, the other including all pairwise interactions with Time. For all three species, both interaction models were much more strongly supported than the models with no interactions (ΔAIC c > 12,766; Appendix S1: Table S4), indicating strong seasonal and diel variation in patterns of habitat selection. We therefore split each species' data for both landscape-and home-range-scale analyses into four subsets for further analysis: dry-day, dry-night, wet-day, and wet-night.

In the final stage of analysis, we selected fixed effects by constructing candidate model sets for each combination of species, scale, season, and time of day. Each set comprised 10 models: (i) intercept only, (ii) all eight additive combinations of the three original predictors, and (iii) a model that included a Mound  Woody interaction term along with the associated main effects. We considered only the Mound  Woody interaction based on our a priori expectation that high woody cover in the matrix might alter mound use (mounds have dense tree cover and thus might be selected less strongly in a woody matrix). All models included a random intercept and random slopes for each predictor, grouped by antelope ID. We used AIC c -based model selection to compare the relative fit of models in each of 24 candidate sets (3 species  2 scales  2 seasons  2 times). We based inferences on the single best-fitting model for each response, which had Akaike weight (w i ) > 0.95 in 21 of 24 sets and w i ≥ 0.71 in the rest, indicating that the top model was always highly likely to be the best in the set. We used 95% Wald-type confidence intervals for the coefficients to assess effect strength (Long et al., 2009). We also report the marginal R 2 (variation explained by fixed effects) and conditional R 2 (variation explained by fixed and random effects) for each model as measures of predictive power (Johnson, 2014;Nakagawa & Schielzeth, 2013).

We conducted a manipulative experiment to test the hypothesis that bushbuck select termitaria thickets for their resources rather than for concealment or shade. We focused on bushbuck, the smallest species, because they should have the highest-quality diets and be most sensitive to predation risk, and thus select most strongly for termite mounds. Moreover, bushbucks' small home ranges made them tractable for this experiment. In July 2015, we used hourly location data from seven collared bushbuck to identify two mounds that were consistently and heavily used by each animal over a 4-week period (n = 14 mounds). We then removed all green (edible) foliage from one of the mounds (selected randomly) up to the 1.5-m maximum browsing height of a bushbuck using pruning shears; we left the other mound as an unmanipulated control.

To test whether our manipulation altered concealment cover, we measured visual obstruction on treated and control mounds using a Robel pole (1.5-m tall, with graduated markings every 10 cm; Robel et al., 1970) both before and after forage removal. We recorded the number of markings that were >50% obscured to an observer standing 5 m away from the mound edge while the pole was held either (i) at the top of the mound or (ii) halfway down the slope of the mound. We repeated this process from each cardinal direction (both before and after manipulation for treatment mounds), yielding eight measurements per control mound and 16 per treatment mound. We compared mean concealment cover among control, pretreatment, and posttreatment mounds using ANOVA, deliberately using measurements (n = 168) instead of mounds as the units of analysis to maximize the odds of detecting an effect (because our aim was to avoid altering cover). We did not manipulate vegetation >1.5-m tall and thus assume that shade and microclimate were unaffected by the treatment.

We used the hourly GPS location data to estimate 95% fixed-kernel UDs (Worton, 1989) for each bushbuck in each of four 1-week periods before and after forage removal. Three of the original focal individuals died before the end of the experiment, leaving n = 4 for analysis (these three animals were not killed by predators and apparently died from dehydration, judging from their atypical directed long-range movements toward permanent water immediately before death). We overlaid weekly UDs for each collared animal onto hand-digitized maps of the treatment and control mounds in each individual's home range. We then calculated the proportion of the volume of each weekly UD that overlapped each of these two mounds. We used the "bushbuck-week" as the unit of replication (total n = 32, each based on ~168 hourly locations: 4 bushbuck  8 weeks, four pretreatment and four posttreatment), which we deem sufficient given that bushbuck traversed their home ranges every ~48 h. We used Student's one-sample t-tests to test the null hypothesis that there was no difference in overlap between bushbuck UDs and (i) treatment mounds before versus after the manipulation (n = 16), and (ii) control mounds before versus after forage was removed from treatment mounds (n = 16).

To evaluate the scaling of behavior with body size, we modified the classic allometric equation, y = a  x b , where y is the response, x is body mass in kg, a is the proportionality coefficient (the intercept at unity), and b is the allometric scaling exponent (the slope of the log-linear regression; Lindstedt et al., 1986;Schmidt-Nielsen, 1984). Most studies of allometry focus on morphological or physiological traits that scale positively with size and cannot have negative values (e.g., metabolic rate, cranial volume). In contrast, behavioral metrics can scale positively or negatively (and convexly or concavely) with size, and responses can be negative (e.g., negative coefficients from a resource selection function indicate avoidance, which is explicitly of interest here). For a consistent approach that would allow a diverse family of functions, we added a constant to give y = a  x b + c, thereby enabling negative y at positive x along with nonzero y-intercepts.

We used this equation to evaluate strength of selection for termitaria at the landscape and home-range scales in each season and time of day. We extracted standardized conditional model coefficients for mound selection using the coef function in R, after accounting for any effects of woody cover and lion use, and then unstandardized and exponentiated those coefficients for inclusion in allometric regressions. To probe the robustness of our inferences, we reiterated these analyses by analyzing use (volumetric overlap between 95% fixed-kernel UDs and mounds; Millspaugh et al., 2006) instead of selection at the landscape scale. We used the same equation to analyze the scaling of home-range area (km 2 ) in each season, mean step length, and mean recursion (revisitation) rate. Unlike for habitat selection and home-range size, we did not conduct separate analyses of step length and recursion for different temporal windows because we had no a priori hypotheses about how these responses should vary across seasons or times of day. Step lengths were calculated as Euclidean distances between successive GPS locations (move package; Kranstauber et al., 2013). Recursion, defined as the number of times an individual returned to a previously occupied site, was calculated by placing a circle of species-specific radius around each GPS location and counting the number of steps that crossed that circle (recurse package; Bracis et al., 2018). We used radii of 150, 425, and 825 m for bushbuck, nyala, and kudu, respectively, following the methods recommended by Fauchald and Tveraa (2003) to determine the scale of area-restricted search behavior for a species (Bracis et al., 2018). Because these species-specific radii inherently account for some of the allometry in space use, we also repeated this analysis using a fixed 150-m radius for each species.

For each allometric regression, we fit nonlinear least-squares models (nls function in R) using a fixed value of a and starting values of b and c that were derived from a preliminary model fit using the Golub-Pereyra algorithm (Golub & Pereyra, 2003). To account for heterogeneity of variance in the regressions for step length and recursion, we weighted those regressions by the inverse of the variance of each metric for each individual. Due to the highly autocorrelated nature of movement data, we used time-series bootstrapping (n = 1000 replicates) to estimate those variances. Scaling relationships were considered allometric when the 95% confidence interval around b (confint2 function in R) did not include 1 (i.e., a sloped line, indicating an isometric relationship) or 0 (i.e., a horizontal line, indicating no relationship). We used AIC c to test whether each model performed better than a null (intercept-only) model, and all but one did; that exception is noted below. We otherwise used the 95% confidence interval (CI) around b for statistical inference.

The Golub-Pereyra algorithm was unable to estimate model parameters in the regressions of home-range size because of two extreme outliers in the dry-season data. These stemmed from kudu that were collared for <3 months during the dry season (two animals died mid-season and their collars were redeployed on two new individuals, which had anomalously large home-range estimates). Short tracking durations can bias estimates of home-range size, particularly for large animals that move long distances (Fleming et al., 2014). We therefore inspected species-level variograms (Calabrese et al., 2016) to estimate the time required to estimate home-range size accurately (2.5 months for kudu and nyala, 1 month for bushbuck). We excluded individuals with tracking periods below these thresholds from analyses of home-range area (only), including the four aforementioned kudu in the dry season and an additional two kudu and 11 bushbuck in the wet season (reflecting high bushbuck mortality in the late dry season, as noted above in relation to the forage-removal experiment). Thus, sample sizes for home-range analysis were 19 bushbuck, 16 nyala, and 18 kudu (n = 53) in the dry season, and 8 bushbuck, 16 nyala, and 16 kudu (n = 40) in the wet season. Because the largest kudu had an anomalously large wet-season home-range estimate, leading to an absurd model, we also reran that model without the outlier (n = 39).

We used DNA metabarcoding to characterize antelope diets. Detailed methods are given in Appendix S2 and broadly followed protocols that we have used to analyze ungulate diets in Gorongosa (Atkins et al., 2019;Becker et al., 2021;Branco et al., 2019;Guyton et al., 2020;Pansu et al., 2019) and Kenya (Kartzinel & Pringle, 2020;Kartzinel et al., 2015Kartzinel et al., , 2019)). We collected fresh fecal samples (n = 52) from the rectums of immobilized antelopes during collaring in 2014 (n = 29; 12 bushbuck, 5 nyala, 12 kudu) and 2015 (n = 23; 7 bushbuck, 6 nyala, 10 kudu) (Appendix S2). We preprocessed samples the same day by transferring homogenized subsamples into tubes containing lysis/preservation buffer, which we vortexed and froze pending transport to Princeton University for DNA extraction and sequencing on an Illumina HiSeq 2500. Analysis focused on the P6 loop of the chloroplast trnL(UAA) intron, a widely used region for plant metabarcoding (Taberlet et al., 2007). After bioinformatic filtering (details in Appendix S2), we identified food-plant sequences by matching them to both an extensive reference library of DNA from locally collected plant specimens (n = 264 sequences, including most of the locally common species; Pansu et al., 2019) and a global reference library derived from the European Molecular Biology Laboratory database (favoring the local library when a perfect match occurred in both). We removed sequences that did not perfectly match any reference sequence, and we merged four sequence pairs that had been assigned the same species name in the local library (retaining the more abundant of each pair for analysis). The remaining sequences were considered molecular operational taxonomic units (mOTUs). After rarefying to the minimum number of sequence reads per sample (n = 6638), the dataset included 164 dietary mOTUs (Appendix S2).

Our analyses are based on relative read abundance (RRA), the proportional representation of each mOTU per sample (Deagle et al., 2019). To estimate forage selectivity, we compared dietary RRA with the relative availability of those plants in the vegetation surveys (details in Appendix S2). We restricted this analysis to eight plant taxa that (i) were matched uniquely to corresponding taxa in the overstory surveys (six species-level and two genus-or subgenus-level identifications) and (ii) averaged >1% RRA in the diet of at least one antelope species across years. These taxa accounted for the majority of diet in each species but did not include several relatively important food plants that we could not identify conclusively in the field and/or match uniquely with an mOTU (notably Diospyros and Combretum spp.). We estimated relative availability by (i) summing the basal area (cm 2 ) of each plant taxon across plots in each habitat type (mound and matrix); (ii) dividing by total area surveyed in each habitat to obtain habitat-specific densities (cm 2 /m 2 ); (iii) multiplying by the proportional areal coverage of mounds (0.011) and matrix (0.989) to obtain scaled habitat-specific densities; (iv) summing across habitats to obtain total availability; and (v) dividing by the cumulative basal-area density of all plant taxa to obtain a proportion. For two congeneric species that we distinguished in the field but had the same DNA barcode (Acacia syn. Vachellia robusta and sieberiana), we summed the availability estimates. Selectivity was calculated for each mOTU in each diet sample as Jacobs's (1974) D, which ranges from 1 to À1 (positive values indicate selection, negative values indicate avoidance). We calculated 95% CIs using the per-sample distribution of D for each antelope species; CIs that did not overlap 0 were considered evidence of significant selection/avoidance. To test our prediction that smaller bodied antelope species eat more of, and select more strongly for, mound-associated plants, we calculated the proportion of habitat-specific density (step iii above) that occurred on mounds and regressed this metric against mean RRA and D for each antelope species.

To visualize overall differences in dietary dissimilarity among samples, species, and years, we used nonmetric multidimensional scaling (NMDS). We tested for significant differences among species and years using factorial permutational MANOVA (perMANOVA) on RRA-weighted Bray-Curtis dissimilarity values among all pairs of samples (vegan package; Oksanen et al., 2019). Subsequently, to test whether diet composition differed between each pair of species based on samples collected in both years, we used pairwise perMANOVA with Benjamini-Hochberg corrections for multiple comparisons. We used factorial ANOVA with Tukey's honestly significant difference (HSD) post hoc contrasts to test for significant differences in dietary richness among species and years.

To estimate nutritional quality of individual antelope diets, we analyzed digestible-energy (DE) and digestible-protein (DP) contents for 25 plant species. We matched these 25 species to 30 mOTUs that collectively accounted for 85% of mean RRA (interquartile range 81%-98%) across all 52 diets. In most cases, we were able to match dietary mOTUs to field-sampled plants with certainty at the species level, but in other cases the matches were only certain at the genus level (details in Appendix S2). We collected young green leaves and stems (petioles) from >3 different individuals of each plant taxon during the mid-dry season (June-August 2016), to match the dietary data. We pooled those samples, dried them to constant weight at 40 C, ground them in a Wiley Mill with a 1-mm screen, and analyzed them for % neutral detergent fiber, % acid detergent fiber, % lignin, % ash, % crude protein, and gross energy at Dairy One Forage Laboratory (Ithaca, NY). From these data, we estimated DE and DP for each plant using the summative equations of Robbins, Hanley, et al. (1987) and Robbins, Mole, et al. (1987). We then calculated weighted averages of DE and DP for each diet using the proportional representation of each food-plant taxon (RRA) as the weighting factor, following the mock-diet approach that we have used previously in this system (Atkins et al., 2019;Branco et al., 2019;Potter et al., 2022). We included samples in this analysis only when >70% of dietary RRA consisted of plants for which we had DE/DP data (n = 47 of 52); the mean RRA of those plants in the remaining diets was 88.1%. We analyzed scaling of DE and DP with body size using the process described above for movement behavior.

Vegetation structure differed sharply between mounds and the matrix. Median canopy cover was 83% on termitaria but 5% in matrix plots (Figure 3a), and median basal-area density was >13-fold greater on mounds (Figure 3b). Of the 45 overstory plant taxa that we identified to genus or species (1594 of 1948 individuals), four were detected only in the matrix and 19 only on mounds. Overstory species composition differed significantly between mounds and matrix and was more dissimilar among matrix plots than among mounds (Figure 3c). Species density was higher on mounds than in the matrix (mean AE SD species per m 2 = 0.06 AE 0.04 on mounds, 0.01 AE 0.01 in the matrix; t = À12.15, df = 117.35, p < 0.001). Rarefaction to common sampling depth (342 individuals) showed that species richness also trended slightly higher on mounds (49.6, 95% CI 41.0-58.1) than in the matrix (43.0, 95% CI 36.6-49.4), although the overlapping CIs indicate that the greater overstory diversity on mounds arose mainly from the difference in stem density. These results collectively support our prediction (P1a).

Species' habitat affiliation strongly predicted the foliar-nutrient concentrations of woody plants (MANOVA effect test, F 16,55 = 5.95, p < 0.0001), consistent with our prediction (P1b). The six mound-associated species were 47% higher in foliar nitrogen, 40% lower in C:N ratio, 100% higher in potassium, 400% higher in sodium and nickel, and >1000% higher in sulfur than the three matrix-associated species; mound-associated species also trended higher in zinc (Figure 3d-f; Appendix S1: Table S5). Boron was the only measured element that was higher (by 26%) in matrix-associated species; the remainder did not differ significantly by habitat affiliation (Appendix S1: Table S5). We found no effect of local growth environment (mound vs. matrix, high vs. low fire frequency, road identity) on any nutrient individually (ANOVA effect tests, p ≥ 0.09) or all collectively (MANOVA effect tests, p ≥ 0.35). These results suggest that mounds filter plant communities in favor of species with high nutrient requirements, contributing to the distinctive composition of termitaria thickets (Figure 3c), but that individuals of a given plant species do not produce more nutritious foliage when growing on termitaria.

As predicted (P1c), all three antelope species selected habitat near termite mounds. This effect was qualitatively consistent across spatial scales, seasons, and times of day, but was most pronounced at the landscape scale, in the dry season, and during daytime (Figure 4a, Appendix S1: Figure S3A; Tables S6 andS7). The landscape-scale models performed well (marginal R 2 = 0.39-0.86, conditional R 2 = 0.94-0.99) and showed that termitaria always had the strongest effect (largest β) on space use by all three antelope species (Appendix S1: Table S6). All species also selected habitat with higher woody cover (except for kudu during the day in the wet season) and lion utilization. The latter result presumably reflects shared habitat preference, but indicates that lion activity was not a strong deterrent (Figure 4b,c, Appendix S1: Table S6). The interaction between mound proximity and woody cover was always negative and significant at the landscape scale (Appendix S1: Table S6), as per our expectation that antelopes should select less strongly for termitaria thickets in areas with high tree density. S1). Mound-associated species had higher (d) percentage nitrogen, (e) sodium, and (f) potassium. Descriptive statistics and statistical tests for all 16 nutrients are given in Appendix S1: Table S5.

Patterns of habitat selection at the home-range scale were broadly concordant with those at the landscape scale, albeit less temporally consistent and with poorer model fits (marginal R 2 = 0.01-0.22, conditional R 2 = 0.23-0.90; Appendix S1: Table S7). Selection for termite mounds was statistically significant in nine of the 12 species  season  time-of-day combinations (and the strongest predictor in six of them) but was more pronounced during the day than at night for each species in each season (Appendix S1: Figure S3A; Table S7). Selection for areas used by lions was more variable within home ranges than at the landscape scale, but those areas were not avoided (coefficients were never significantly negative; Appendix S1: Figure S3B), bolstering the inference that lion activity (and the low level of lion predation on Tragelaphus spp.; Bouley et al., 2021) was not a major source of perceived risk. As at the landscape scale, all three species selected for high woody cover within their home ranges, except for kudu during the day in the wet season (when they significantly avoided woody cover; Appendix S1: Figure S3C). The interaction between mound proximity and woody cover was significantly negative in seven of the 10 models in which it was included, again broadly mirroring the pattern at the landscape scale (Appendix S1: Table S7).

Removing available forage from termite mounds did not significantly reduce concealment cover (Figure 5a) but did reduce use of mounds by bushbuck in the 4 weeks following forage removal, whereas use of control mounds did not change (Figure 5b,c). These experimental results, consistent with our prediction (P1c), show that forage is a strong driver of habitat selection by bushbuck in Gorongosa and a sufficient explanation for their use of termitaria, irrespective of any contributing effects of concealment cover or shade.

As predicted (P2a), strength of selection for termitaria generally declined nonlinearly with body mass (strongest in bushbuck and weakest in kudu) albeit with considerable scatter arising from individual-level variability, especially in bushbuck (Figure 6; Appendix S1: Figure S4; Tables S8 andS9). The allometry of mound selection was most pronounced in the dry season and during the day, especially at the home-range scale, where there was no real relationship during the wet season (the scaling exponent was indistinguishable from 0 in both wet-season b) Bushbuck and nyala selected for woody cover in all intervals, whereas kudu selected for woody cover only during the dry season. (c) All species selected significantly for areas that were used by lions, which probably reflects shared habitat preferences rather than causal linkage, but indicates that antelopes were not deterred from using these areas. Analogous results for the home-range scale are given in Appendix S1: Figure S3.

models at the home-range scale; Appendix S1: Table S9). The shape of these regressions was consistently concave-up, decreasing in the dry season but mixed in the wet season; note, however, that wet-season analyses were compromised by the small sample of bushbuck (n = 8, vs. 19 in the dry season). Analyzing mound use instead of mound selection at the landscape scale yielded very similar results (Appendix S1: Figure S5). We found no coherent support for our prediction (P2b) that large-bodied antelopes would exhibit weaker selection for mounds in the wet season due to increased forage availability in the matrix. At the landscape scale, all species exhibited weaker selection for mounds in the wet season, which is consistent with our intuition that termitaria should be less attractive in the wet season, but not with the predicted size structure of that effect; although strength of selection for mounds did decline more steeply at high body-mass values in the wet season than the dry season (Figure 6), seasonal differences were generally nonsignificant at the species level and were, if anything, greater in smaller bodied species (bushbuck and nyala; Figure 4a). At the home-range scale, only bushbuck exhibited any decrease in selection for mounds from dry to wet season, opposite to our prediction (Appendix S1: Figures S3A andS4). Home-range sizes were larger in the wet season and increased as a concave-up function of body mass in both seasons (Figure 7a; Appendix S1: Figure S6; Table S10). In the wet season, our fitting process yielded a ridiculous model of home-range size (coefficient a = 10 À44 and scaling exponent b = 21.63, roughly a right angle), which arose from an anomalously large home-range estimate for the largest kudu; excluding this outlier resulted in a more plausible model with b that remained both significantly >1 and significantly higher than that in the dry season (Appendix S1: Figure S6; Table S10). Step length scaled positively and nonlinearly (and concave-down) with size (Figure 7b), while recursion scaled negatively and nonlinearly (and concave-down; Figure 7c); re-running the recursion analysis with a common radius for all three species (150 m) amplified the differences among species and resulted in an even more pronounced allometric relationship. Areas near termite mounds were revisited most frequently by bushbuck, less so by nyala, and least so by kudu (Figure 7d). These results accord with our prediction (P2c) and indicate a greater tendency of small-bodied antelopes to concentrate foraging in smaller areas, especially around termitaria.

Body mass (kg) y = -0.33x 0.13 + 1.67 y = -1.51× 10 -8 x 3.08 + 1.08 y = -1.02 × 10 -6 x 2.23 + 1.08 y = -6.77x 0.0049 + 8.00 S8; the 95% CIs around b never overlapped 1 or 0, indicating significant nonlinearity in all temporal windows. Analogous results for the home-range scale are given in Appendix S1: Figure S4.

Body size, diet composition, and dietary niche (Hypothesis 3) S10; regressions for step length and recursion were weighted by the inverse of the variance for each individual to account for autocorrelation. (d) Number of revisitations (10,000 randomly chosen points per species) as a function of distance from the site centroid to the nearest termite mound; x-axes are truncated at 400 m (tails extend to 3000 m). Smaller bodied species more frequently revisited sites centered on or near mounds.

for kudu). In contrast, Combretum spp. were scarce on mounds-the matrix supported seven-fold higher basal-area density across all taxa and 82%-100% of individuals for the two most abundant species (Appendix S1: Table S1)-and were eaten in quantity only by kudu (14% mean and up to 83% individual RRA, vs. 2% mean and up to 13% individual RRA for bushbuck and nyala).

The eight plants for which we calculated selectivity did not include Diospyros or Combretum spp. but nonetheless accounted for 65.8%, 50.8%, and 60.4% of mean RRA for bushbuck, nyala, and kudu, respectively. All of these plants occurred at higher densities on mounds than in matrix plots (Figure 8a), although the proportion of their total abundance that occurred on mounds ranged from <2% to 100% (Figure 8b); this discrepancy reflects the low coverage of termitaria in the landscape (1.1%), such that availability of a species can be higher in the matrix even when mounds support denser concentrations. Bushbuck selected three of these taxa relative to their availability, nyala selected four, and kudu selected two (Figure 8c). The top food plant for both bushbuck and nyala (B. discolor) occurred exclusively on mounds, whereas the least mound-associated taxon (Acacia syn. Vachellia robusta/sieberiana) was the only one avoided by all three species (Figure 8c). Mean RRA and selectivity for these taxa increased significantly with the proportion of their total abundance that occurred on mounds for bushbuck (RRA: R 2 = 0.62, F 1,6 = 9.64, p = 0.021; selectivity: R 2 = 0.68, F 1,4 = 8.33, p = 0.045) and nyala (RRA: R 2 = 0.64, F 1,6 = 10.79, p = 0.017; selectivity: R 2 = 0.59, F 1,6 = 8.81, p = 0.025) but not for kudu (RRA: R 2 = 0.09, F 1,6 = 0.62, p = 0.46; selectivity: R 2 = 0.15, F 1,4 = 0.70, p = 0.45). These results support our prediction (P3a) that mound-associated woody plants are prevalent in each species' diet but account for a greater share of diet in (and are selected more strongly by) smaller bodied species.

As predicted (P3b), diet composition differed significantly among species and these differences were aligned with body size, with bushbuck and kudu being most dissimilar (Figure 9a). Post hoc contrasts of dietary dissimilarity were statistically significant between each pair of species after correcting for multiple comparisons S11. DP and DE also differed significantly among species (ANOVA: DP F 2,43 = 11.32, p = 0.0001; DE, F 2,44 = 17.70, p < 0.0001); lowercase letters in the legend keys indicate statistically significant differences between species in Tukey's HSD post hoc tests (all p < 0.002). One outlying kudu (KD_13A) with extremely low DP (4.83%) was omitted from that analysis, mainly for graphical convenience; exclusion of this outlier did not alter the significance of the species-wise contrasts or the functional form of the regression (Appendix S1: Table S11).

(bushbuck vs. kudu pseudo-F 1,39 = 20.74, R 2 = 0.35, p ≤ 0.0015; nyala vs. bushbuck pseudo-F 1,28 = 2.12, R 2 = 0.07, p = 0.0280; nyala vs. kudu pseudo-F 1,31 = 13.82, R 2 = 0.31, p ≤ 0.0015). Dietary richness was significantly higher in kudu and nyala (least-squares means 61 and 60 mOTUs, respectively) than in bushbuck (51 mOTUs) and was higher in 2015, the drier year (whole-model F 5,46 = 12.88, p < 0.0001; species F 2,46 = 11.44, p < 0.0001; year F 1,46 = 22.96, p < 0.0001; species  year F 2,46 = 3.09 p = 0.06; Tukey's HSD contrasts for bushbuck p ≤ 0.004). Dietary differentiation arose from divergent patterns of selectivity (Figure 8) in addition to the aforementioned differences in consumption of Diospyros and Combretum spp.

For example, the top kudu food (C. kirkii, 32.9% RRA) was significantly avoided by bushbuck and nyala (<7% RRA); the top food for bushbuck and nyala (B. discolor, >20% RRA) was essentially uneaten by kudu (0.6% RRA); and the second-ranked and strongly selected bushbuck food (T. capitata, 15.2% RRA) was only marginally selected by nyala (6.9% RRA) and was avoided by kudu (1.2% RRA). We found mixed support for our prediction (P3b) that diet quality would scale negatively with body size (Appendix S1: Table S11). As expected, DP decreased allometrically with body mass and was significantly higher in bushbuck than nyala and kudu (Figure 9b). Unexpectedly, DE increased allometrically with body mass and was significantly higher in kudu than bushbuck and nyala (Figure 9c). Despite the contrasting direction of these relationships, they contributed to differentiating the nutritional niches of each species: bushbuck diets had high DP but low DE, nyala had low DP and low DE, and kudu had low DP but high DE (Figure 9b,c).

We found strong support for seven of our nine hypothesis-driven predictions, mixed support for one, and no support for another (Table 1). Altogether, our results show that foraging behaviors scale allometrically with body size in ways that collectively differentiate the use of space and food by sympatric, congeneric ungulates. Spatial and dietary niche differentiation stabilizes coexistence in animal communities (Pringle et al., 2019) and is often associated with body-size differences among closely related species-notably, for example, in the adaptive radiations of Darwin's finches (Grant & Grant, 2008), sticklebacks (Schluter, 2000), and Anolis lizards (Losos, 2009)-suggesting that selection on body size is important in both the origin and maintenance of species diversity. Our results are consistent with this view and show that allometric scaling of resource selection is a behavioral mechanism through which size-structured niche differences emerge, thereby linking the allometry of behavior (Dial et al., 2008) with species coexistence. Theoretical models of body-size-mediated coexistence have incorporated assumptions about the allometry of nutritional requirements and space use in patchy landscapes (Bagchi & Ritchie, 2012;Basset & Angelis, 2007;Hopcraft et al., 2010;Ritchie & Olff, 1999), but few studies have empirically tracked how these relationships lead to resource partitioning among free-ranging animals.

Termitaria supported distinct assemblages of nutrient-rich plant species, and all antelope species selected for mound-associated space and food. Strong effects of fungus-farming termites on plant species composition, productivity, and nutrient content-arising from alteration of soil structure, chemistry, and moisture-appear to be nearly universal in African savannas and tend to decay nonlinearly with distance from mound centers (Baker et al., 2020;Joseph et al., 2013;Moe et al., 2009;Pringle et al., 2010;Seymour et al., 2014;Sileshi et al., 2010;Sileshi & Arshad, 2012). Many studies have also found that mound-associated plants are heavily used and affected by herbivores (Grant & Scholes, 2006;Holdo & McDowell, 2004;Joseph et al., 2018;Levick et al., 2010;Loveridge & Moe, 2004;Mobaek et al., 2005;Okullo et al., 2013;Seymour et al., 2016), but this effect is not universal and several studies have reported the opposite pattern (Davies, Levick, et al., 2016;Muvengwi et al., 2013Muvengwi et al., , 2019;;O'Connor, 2013;Van der Plas et al., 2013).

Although our tests of Hypothesis 1 are largely consistent with the literature on termite-herbivore interactions in savannas, they also differ from previous work in several important ways. First, prior studies have assessed mound use at small scales (from individual mounds up to a few km 2 ) by counting dung, observing animals, or quantifying herbivore impacts on plants. Our use of GPS movement data together with high-resolution remotely sensed imagery enabled us to quantify strength of selection for mounds at both large spatiotemporal extents and fine spatiotemporal grains while controlling for other factors that shape habitat use. This, in turn, equipped us to assess context-dependency in the strength of selection for mounds (typically stronger at the landscape scale, in the dry season, during the day, and in areas with lower woody cover in the matrix) and to conclude that mound proximity remained the overriding determinant of habitat selection even after accounting for effects of woody cover and predator activity.

Second, ours is the first study to experimentally test the mechanistic basis of mound use by large herbivores. The role of forage availability is often assumed, but other factors could also promote the frequently observed pattern of higher foraging intensity on mounds. For example, mounds might limit predation risk by providing elevated vantage points for vigilance and/or concealment for "hider" species such as bushbuck (Anderson et al., 2016;Atkins et al., 2019). Mound thickets also have cool microclimates and might thus be thermal refuges (Joseph et al., 2016). Our data show that forage availability alone was sufficient to drive patterns of mound use, irrespective of risk or shade. We were only able to test this hypothesis for bushbuck (due to their small home ranges and the effort required to remove forage from mounds), but several lines of evidence suggest that selection for mounds by nyala and kudu is likewise driven by resource availability: (i) all three species disproportionately ate plants that occurred at high densities on mounds; (ii) a primary role of microclimate would suggest stronger selection for mounds in the hot wet season, when mammals are most likely to overheat, whereas we found the opposite; and (iii) any risk-driven effects should be strongest for bushbuck, which rely on crypsis (Estes, 2012) and are vulnerable due to their small size and solitary habit (Sinclair et al., 2003).

A caveat to this conclusion is that predation on these species was very low during our study. This low-risk environment was key to our study design because it enabled us to isolate allometries of resource-acquisition behavior, but it may also have led to patterns of habitat use that differ from those in more intact systems. Indeed, we have found that bushbuck in Gorongosa relaxed their risk sensitivity and expanded into open habitat following the extirpation of leopard and wild dog (Atkins et al., 2019). It is therefore possible that termitaria are more important for risk avoidance in other systems and may be increasingly important in Gorongosa as predation pressure intensifies. Wild dogs were reintroduced in 2018 and have fed heavily on bushbuck (Bouley et al., 2021), and leopards have been introduced since that time. If termitaria provide safe spaces for small-bodied species in particular, then predators might accentuate the observed allometries and niche differences. Alternatively, if mounds confer safety for antelopes of all sizes, then risk might collapse niche differences as species converge on the same refuges (Pringle et al., 2019). Our study provides a baseline that future work can use to test these possibilities as carnivores recover in Gorongosa. However, anecdotal observations in 2021-wild dogs hunting around mounds to flush bushbuck, fewer bushbuck atop mounds than in previous years-suggest that mounds may not be particularly safe after all, and that increased predation pressure might actually decrease selection for mounds. The woody cover on termitaria is also a double-edged sword in that it provides concealment for predators as well as for prey (Davies et al., 2016b;Ford et al., 2014;Valeix et al., 2009), a likely explanation for our finding that antelopes and lions exhibited similar patterns of habitat use.

Ultimately, our results show clearly that food availability was a major determinant of mound use (and habitat selection more broadly) by our focal species, and previous work suggests that this is a common pattern across many savanna herbivore species and mound types. Still, risk avoidance and favorable microclimate are not mutually exclusive with high resource availability, and all three factors may contribute to explaining mound use by savanna ungulates. Future work aimed at parsing these potentially complementary mechanisms would enrich our understanding of the ways in which termites shape patterns of biodiversity, animal behavior, and ecosystem function in savannas (Castillo Vardaro et al., 2021;Pringle et al., 2010).

All three antelope species exhibited broad similarities in habitat use but diverged in strength of selection for habitat features and other aspects of movement behavior. In general, these metrics scaled nonlinearly with body mass across species. However, these relationships were rarely discernible within species (where individual variability often matched the range of interspecific variation), and allometric regressions took all possible functional forms (increasing, decreasing, concave-up, concave-down). Thus, our study strongly supports the proposition that movement behavior scales allometrically with body size, but explaining the full suite of relationships in terms of underlying morphological and physiological traits will require further research.

Strength of selection for termitaria decreased with body size in all spatiotemporal windows except for at the home-range scale during the wet season at night. Overall, selection for mounds was stronger, and allometric scaling more pronounced and consistently concave-up, during the dry season. This pattern accords with the logic underlying Hypothesis 2b, that selection for mounds is stronger when food availability is limited in the matrix, yet we found no support for the expected size-dependence of this effect (although our ability to detect it was weakened by the limited sample of bushbuck in the wet season). To our knowledge, only one study has explicitly tested the possibility that smaller herbivore species use termite mounds more intensively. Mobaek et al. (2005) found no relationship between body size and mean foraging distance from mounds among five species (oribi, bushbuck, impala, topi, and warthog) in Uganda. However, that study did not control for relatedness or even broad dietary differences (topi and warthog are grazers, oribi and impala are mixed feeders, bushbuck are browsers). We focused on congeneric browsers to minimize any such confounding effects of phylogeny and associated trait differences. Even our carefully chosen set of focal species exhibit trait differences that might add noise to allometric signals in behavior. For example, bushbuck are solitary, whereas nyala and kudu forage in small groups. This difference in sociality, unaccounted for in our analyses, could produce variation in group-level resource requirements and associated patterns of space use. That we still detected clear allometric signals in selection for termitaria suggests that any such effects are not so strong as to override the influence of body size.

Intraspecific variation in strength of selection for mounds was greatest among bushbuck, at least during the dry season (individuals of all species varied widely in their generally weaker selection for mounds during the wet season). The individual variability among bushbuck may arise in part from their small home ranges (dry season mean = 2.7 km 2 ) coupled with variation in the spatial distribution of mounds, which are locally overdispersed but regionally patchy (Figure 1c). As a result, some bushbuck simply might not have many mounds to select. Indeed, the density of mounds within bushbuck home ranges ranged widely (dry season range 66-231 km À2 ) and was more variable than for other species (Appendix S1: Figure S7). Moreover, bushbuck and nyala with lower mound densities in their home ranges generally exhibited weaker selection for mounds (Appendix S1: Figures S8 andS9). Such functional responses in selection (changes in selection as availability changes; Godvik et al., 2009;Mauritzen et al., 2003) may represent an adaptive strategy for coping with resource heterogeneity across space and time, but are most often negative (i.e., selection for a favorable habitat increases when that habitat is rare; Anderson et al., 2012). The positive functional response to variation in mound density is thus unusual and may reflect the existence of flexible behavioral strategies depending on local ecological context: when termitaria are too scarce to sustain a "mound-hopping" strategy, individuals may shift to a foraging strategy that enables more efficient use of matrix-associated plants. This interpretation resonates with our previous finding that, in the absence of strong predation risk, many Gorongosa bushbuck have expanded into open floodplain habitat with few termite mounds but with otherwise abundant high-quality forage (Atkins et al., 2019). Some of the residual variation observed in all species might also stem from state-dependent behaviors by individuals that differed in nutritional condition, lactation, and/or gestation status (Owen-Smith et al., 2010). Resolving individual-level drivers of behavioral variation may help to reduce noise at the species level and make allometric signatures in behavior easier to detect.

Our results for home-range area, step length, and recursion also support the hypothesis that the concentration of food on termitaria leads to allometric scaling of foraging intensity, with small animals using smaller areas more vigorously. The relationship between home-range area and body size in mammals is the focus of a large body of literature and a long-standing debate over whether and when scaling should be sublinear, linear (isometric), or superlinear (Calder III, 1983;Harestad & Bunnell, 1979;Haskell et al., 2002;Jetz et al., 2004;Kelt & Van Vuren, 2001;Mcnab, 1963;Tamburello et al., 2015). Our finding that home-range size increased as a concave-up function of body mass accords with the most recent and taxonomically inclusive analysis, involving 61 globally distributed mammal species (Noonan et al., 2020); indeed, our scaling exponent in the dry season, when sample size and tracking duration were greatest (b = 1.40), almost perfectly matches the one reported for herbivores in that study (b = 1.38). This consistency is noteworthy in light of the observation that ecological patterns in global syntheses may not "scale down" to the community level (Messier, Lechowicz, et al., 2017;Messier, McGill, et al., 2017). Although our analysis of home-range scaling in the wet season is less robust due to the smaller number of individuals, shorter tracking durations, and several outlying estimates, the functional form of the relationship was consistent.

We likewise found that large body size was associated with larger step lengths (sublinearly, unlike home-range size), lower rates of site revisitation (recursion), and less frequent revisitation of sites near termite mounds. Together with our results for home ranges, these patterns indicate higher levels of area-restricted search behavior in smaller bodied species. Smaller step lengths in bushbuck relative to nyala and kudu may reflect the mechanical and metabolic constraints on locomotion imposed by small body size (Cloyed et al., 2021;Peters, 1983;Portalier et al., 2019). However, bushbuck revisited resource hotspots more frequently than the larger bodied species, suggesting a potential role of spatial memory in the allometric scaling of foraging behavior. Theoretical models suggest that memory-based movement could explain the higher frequencies of spatial segregation-often assumed to reflect territoriality-observed in small-bodied antelope species (Riotte-Lambert et al., 2015). This possibility is intriguing in relation to bushbuck, which are solitary but not highly territorial and which, according to Estes (2012, p. 173) "should not be considered antisocial but rather as loosely and casually sociable." Empirical research on spatial memory in ungulate foraging remains limited and is an exciting avenue for future work.

Body size, diet composition, and dietary niche (Hypothesis 3)

As with movement behavior, superficial similarities in resource use among our three sympatric species concealed significant dietary niche differences, as indexed by the plant taxa they selected and the dissimilar composition and nutritional quality of their diets. Multiple lines of evidence suggest that allometric scaling of movement behavior (and selection for termitaria especially) is directly connected to these dietary differences, consistent with work showing that interspecific differences in space use and movement often correspond with differences in diet (Fleming, 1991;Sailer et al., 1985). The plant that accounted for a third of bushbuck diets and a fifth of nyala diets, B. discolor, occurred only on mounds and was rare in the landscape. By contrast, the plant that accounted for a third of kudu diets, C. kirkii, was abundant in both habitats and was one of the dominant woody plants in the landscape (150-fold higher density than B. discolor). Bushbuck and nyala selected more strongly for plants that were more strongly mound affiliated, whereas kudu did not. Moreover, only kudu ate substantial amounts of Combretum spp., which were unique among the major food plants in occurring almost exclusively in the matrix.

The large literature on resource partitioning in African savanna ungulates has emphasized two axes of differentiation. One is the proportion of grass versus browse in the diet. By this criterion alone, our three species were equivalent: all were strict browsers with diets comprising <1% grass on average (range 0%-3% and median 0.3% across all 52 individuals). Studies elsewhere have found that nyala eat grass more often than bushbuck or kudu (Codron et al., 2007), but we found no such signal. Only by characterizing the taxonomic makeup of diets were we able to detect interspecific differences in diet composition, richness, and selectivity, which underscores the utility of DNA metabarcoding for resolving cryptic niche differences (Pringle & Hutchinson, 2020). These results for Tragelaphus spp. in Mozambique are analogous to our previous finding that sympatric zebras (Equus spp.) in Kenya partitioned plant taxa despite both occupying the same position on the grazer-browser spectrum (Kartzinel et al., 2015).

The other canonical axis of diet differentiation in ungulates is forage quality, and there, body size and allometric scaling are thought to play a crucial role. The Jarman-Bell hypothesis holds that large-bodied species eat lower quality diets than small-bodied ones (Bell, 1971;Geist, 1974;Jarman, 1974). Allometric scaling of basal metabolic rate (≈ M 0.75 ) and food-retention time (thought to be ≈ M 0.25 ) was central to this idea as first conceived and later developed (Demment & Van Soest, 1985;Illius & Gordon, 1987): smaller herbivores require higher quality diets because of their high mass-specific metabolic rates, whereas larger herbivores can tolerate lower quality food by retaining it longer and digesting it more fully. In the last decade, however, evidence has mounted that, although inverse relationships between size and diet quality are empirically common, they do not emerge straightforwardly from physiological allometries (Clauss et al., 2013;Müller et al., 2013;Potter & Pringle, 2022;Steuer et al., 2014). In particular, there is little evidence that larger herbivores have greater digestive ability (retention time ≈ M 0 ; Steuer et al., 2011). Rather, larger herbivores might have lower quality diets on average simply because they eat more food (intake ≈ M 0.9 in ruminants; Hackmann & Spain, 2010) and are less able to efficiently select or maintain themselves only on high-quality forage, which tends to be sparse. Clauss et al. (2013, p. 13) concluded "that animals of any size can use diets of any quality" and that nutritious diets should be desirable for herbivores of all sizes, but that "larger herbivores are (mostly) confined to low-quality diets," whereas small ones can afford to be picky and generally will be, except when "excluded from a certain range of plants or plant parts because of physical limitations."

Against this backdrop, our finding that DP scaled negatively with size (higher in bushbuck than nyala and kudu, as expected), whereas DE scaled positively (higher in kudu than bushbuck and nyala, contra expectation) is intriguing, especially given that DE and DP were positively correlated across plant taxa (R 2 = 0.19, F 1,21 = 4.85, p = 0.039). This result shows that empirical tests of the Jarman-Bell hypothesis are sensitive to which aspect of "quality" has been assessed. Previous work has used diverse proxies for diet quality, including dietary crude protein, fecal N, consumption of nonstem material, short-chain fatty acid concentration, fermentation rate, and preference for burned patches, and these metrics showed broadly concordant negative allometric scaling with size (Clauss et al., 2013). Yet, few prior studies have analyzed multiple metrics of quality simultaneously for multiple free-ranging sympatric species. Seasonal dynamics of fat and protein reserves play key roles in the fitness and life histories of long-lived mammals (Monteith et al., 2013). Accordingly, both DP (which influences lean muscle mass and associated protein reserves; Monteith et al., 2013) and DE (which affects fat accretion; Cook et al., 2004) are thought to be useful metrics of ungulate diet quality. We hypothesize, however, that DP is a particularly good metric in the tropics, where ungulates often rely on an income-breeding strategy (Jönsson, 1997;Sainmont et al., 2014;Stephens et al., 2014) and do not accumulate large fat reserves even when conditions are favorable.

The observed differentiation in diet quality can be traced to differences in the use of plant species that differed in their affinity for termitaria. The top kudu food (C. kirkii, widespread on and off mounds) had the highest DE of any major food plant measured but lower-than-average DP, while the third-ranked kudu food taxon (Combretum spp., matrix affiliated) had above-average DE but extremely low DP, helping to explain why kudu diets were high in DE but low in DP. Bushbuck ate high proportions of B. discolor, a protein-rich but energetically mediocre plant found only on mounds, helping to explain their high DP. Nyala also ate substantial amounts of B. discolor, but that proportion was lower than for bushbuck and offset by a higher proportion of Diospyros spp., which had very low DP and middling DE. Although these results are consistent with the idea that body size does not physiologically necessitate diets of a particular quality, we can only speculate about why DE and DP scaled in opposite directions with body mass. It might be, for example, that mound-restricted plants such as B. discolor are insufficiently abundant to sustain the higher absolute requirements of kudu (especially in competition with bushbuck and nyala), but that kudu compensate by ranging farther and eating comparatively energy-rich plants (C. kirkii, Combretum spp.) that are more widely distributed in the landscape. Regardless of its mechanistic basis, the allometric scaling of diet quality served to differentiate each species' nutritional niche, even if not in the expected fashion: nyala, the medium-sized species, had the lowest-quality diet when considering DP and DE together.

We acknowledge that our approach to estimating diet quality, like all methods available for wild ungulates, is imperfect and subject to several caveats. The Jarman-Bell hypothesis is generally interpreted with respect to the consumption of plant tissues that differ in quality, such as fruits and young shoots versus old leaves and stems (du Toit & Olff, 2014;Hopcraft et al., 2010). DNA metabarcoding has the crucial advantage of high taxonomic resolution, but it cannot reveal which tissues were eaten, nor at what height; our inferences are based on differences in the plant species eaten by each antelope and the (substantial) interspecific variation in the quality of the green leaves and stems that we sampled. This caveat, also inherent in some of our prior work (Atkins et al., 2019;Becker et al., 2021;Branco et al., 2019), cannot qualitatively affect our core conclusion of size-structured dietary niche differences: interspecific variation in the use of particular tissues or browsing heights within a plant would only augment the differences already evident in taxonomic diet composition, so our approach is conservative in this regard. Similarly, any such differences should only accentuate the negative scaling of DP with size, as small-bodied species are more able to consistently select high-quality tissues. It is possible that our approach underestimates DE for small, relative to large, species; for example, bushbuck will eat fruits when available (Estes, 2012;Kingdon, 2015) and several termitaria species (e.g., B. discolor) have palatable and presumably energy-rich fruits. However, we sampled diets in the mid-dry season (June-August), which is not a peak flowering, fruiting, or leafing time in this system (Tinley, 1977).

Another potential source of error in the DE/DP analyses is imprecision and gaps in our matching of dietary mOTUs with plants sampled in the field (Appendix S2), which in several cases assumed that congeneric plants have similar DE/DP (which is not always true). While we were able to match most major forages with a high degree of certainty and think it unlikely that the remaining gaps would qualitatively affect our inferences, these gaps nonetheless highlight the importance of local reference libraries for maximizing the value of metabarcoding data. Our local library for Gorongosa plants is extensive, but not comprehensive (Pansu et al., 2019), and future additions could further refine the taxonomic assignment of sequences obtained in this study. Such refinements will not resolve cases where closely related plants share the same barcode but, in some cases, mOTUs were only identified to genus or family because the corresponding plant species was either not collected or not identified beyond genus in the reference library (Appendix S2). Future studies might also include field observations of foraging animals to gain complementary data on differences in the use of plant parts (leaves, stems, fruits) and browsing height, which may enable herbivores to partition food resources to an even greater degree than revealed by taxonomic composition alone (du Toit, 1990). Our diet samples were collected in the mid-dry seasons of successive years, when we expected moderate resource partitioning (stronger than the wet season, weaker than the late dry); sampling across seasons, as we did for movement, would illuminate the extent to which diets converge and diverge depending on resource availability.

We emphasize one caveat to fitting allometric models using nonlinear regression. This approach was needed to accommodate a variety of scaling relationships and, in nearly all cases, the models fit trends that were visually obvious in the data. Importantly however, this method precludes the comparison of our scaling exponents (b) with those from the standard allometric-regression model, y = a  x b , except in the few cases (e.g., home-range size) in which our model converged to something very similar to the standard model. Our addition of the constant c to these models means that we can obtain similar functional forms in different ways depending on the signs of a and b. For example, a concave-down increasing function, as found in many established allometric relationships, occurs when 0 < b < 1 in the standard equation (and in ours provided a > 0), but can also occur in our models when both a and b are negative (in which case c enables positive values of y). We therefore consider our models to be valid descriptions of the functional form of the relationship between body mass and each response variable, but we do not attribute inherent meaning to the fitted values of these terms and caution against extrapolating the models beyond the range of body masses in our data.

The variable functional forms of our models also make it hard to propose unifying mechanistic explanations based on underlying morphophysiological traits, and indeed suggest the interplay of multiple size-linked traits. The prospect of being able to explain ecological phenomena in terms of biophysical first principles is alluring (Brown et al., 2004), but beyond our reach here. We hope that future work succeeds in discovering size-structured traits and processes that underpin our results. We also note that, even for the best-studied relationships in our study such as those between body size and home-range area (Haskell et al., 2002;Jetz et al., 2004;Noonan et al., 2020) and diet quality (Clauss et al., 2013;Müller et al., 2013), there is no firm consensus about the exact functional form of the relationships, much less their mechanistic basis.

We conclude that allometric scaling of behavior leads to size-structured niche differences in a set of closely related sympatric ungulates. These differences are not superficially obvious, given the broad interspecific similarities in space and resource use, but were detectable via high-resolution analysis of space use and diet. Niche differences mitigate competition and are necessary for coexistence but not sufficient, because large discrepancies in competitive ability can overwhelm stabilizing niche differences (Broekman et al., 2019;Chesson, 2000). Our study does not establish that the observed niche differences enable coexistence of these three species, nor even that they are stably coexisting (Chesson, 2000). However, their historical co-occurrence in Gorongosa and elsewhere (Tello & Van Gelder, 1975;Tinley, 1977) coupled with their similar postwar recovery rates suggest that these species can coexist in Gorongosa, and we hypothesize that the niche differences documented here enable that outcome. It is tempting to speculate that allometries of behavior (Dial et al., 2008) underpin coexistence of close relatives worldwide, where size-based niche differences and character displacement/release are common (Bagchi & Ritchie, 2012;Basset & Angelis, 2007;Ritchie & Olff, 1999).

Our findings have regional implications for understanding the dynamics of spiral-horned antelopes. Competition between bushbuck and nyala in particular is a long-standing management concern in KwaZulu-Natal, South Africa, where bushbuck have declined to the point of extirpation in several places in concert with increasing nyala density. It is widely perceived that these declines stem from interspecific competition (Coates & Downs, 2005;Fay & Greeff, 1999). One study concluded that nyala had competitively excluded bushbuck in the Ndumo Game Reserve (Ehlers Smith et al., 2020); another found this explanation insufficient for bushbuck decline in Hluhluwe-iMfolozi Park (Owen-Smith et al., 2017), based in part on the theoretical expectation that smaller ungulates should outcompete larger ones (Bagchi & Ritchie, 2012;du Toit & Olff, 2014). Our study affirms overlap in resource use between nyala and bushbuck-we found no evidence for exclusive resources-but also shows that they partitioned shared resources in hitherto undetected ways. Thus, in Gorongosa, stabilizing niche differences were larger, and the degree of competitive equivalence needed to permit coexistence smaller, than expected from prior work. If all Tragelaphus spp. eventually reach stable densities in Gorongosa after having all increased from low abundance, then we would view that as compelling evidence that these niche differences do indeed enable stable coexistence. Definitively linking niche differences to competition and coexistence outcomes remains a crucial empirical challenge for ungulates (Prins, 2016) and animals in general (Ritchie, 2002), which are less tractable than plants for the experiments used to test coexistence theory (Broekman et al., 2019;Grainger et al., 2019).

A final conclusion of our study is that fungus-farming termites play a vital role in governing savanna dynamics by creating spatial heterogeneity that serves as a template for niche differentiation in plants and animals. The published literature on termites as ecosystem engineers in savannas attests to their role in sustaining biodiversity and ecosystem functions (Castillo Vardaro et al., 2021;Joseph et al., 2011Joseph et al., , 2013;;Pringle et al., 2010). That they may also help to anchor the two defining attributes of African savannas-tree-grass coexistence (please refer to Appendix S1: Figure S1) and ungulate coexistence-is an intriguing prospect that warrants further investigation.