Out of scale out of place: Black rhino forage preference across the hierarchical organization of the savanna ecosystem

The successful conservation plans of megaherbivores necessitate precisely characterizing their ecological needs in order to optimize reproduction rates and reintroduction plans. The black rhino (Diceros bicornis L.) is among the most endangered species of megaherbivores in Africa and its conservation relies on nature reserves that are bound and habitat‐restricted. Therefore, identifying the optimal amount of space this species needs and the factors driving its habitat use are crucial for establishing reserve priority plans. Knowing that forage selection is an important component linked to herbivore spatial distribution, we combined 5 years of sightings data with observations of rhinos' vegetation type and forage preferences to address their forage selection across multiple spatial scales. We found that black rhinos' spatial distribution was negatively associated with ecosystem productivity, but positively associated with specific vegetation types that contain highly preferred, chemically distinct, plant species. Black rhinos thus occupy their habitat across space and time through selective foraging on preferred plants.

the worldwide population since 1960 has declined by an estimated 98% (Emslie & Adcock, 2016). Poaching for rhino horn is the species' biggest threat, but not the only one. The lack of habitat availability also threatens their survival (Emslie & Adcock, 2016). At present, wild black rhino populations are limited within bounds of protected areas, which are often enclosed and thus restrict migration (Landman et al., 2013). In this context, the number of individuals within protected areas needs to be managed to remain just under carrying capacity in order to maintain and increase the global black rhino population. This strategy comprises selecting animals from populations which have a positive growth rate and translocating them to alternative areas that have potential to sustain new populations (Linklater et al., 2012).
Apart from security concerns, the main challenge for successful translocation plans, resides in finding habitats that have all necessary components for maximizing the species' growth rate (Balfour et al., 2019;Odendaal-Holmes, Marshal, Parrini, & Parrini, 2014). In characterizing suitable habitats and subsequently estimating a habitat's potential carrying capacity, many factors need to be considered, including climatic conditions and seasonal variation, plant productivity, forage selectivity, and topography (Adcock, 2001). However, information precisely identifying black rhino habitat use is still generally not available. Black rhinos, like most large mammalian herbivores, forage across several temporal and spatial scales. In heterogeneous landscapes, they select habitat patches mainly to consume specific plants to maximize nutrient and energy intake (Owen-Smith, Fryxell, & Merrill, 2010), however, other factors such as home ranges (Mitchell & Powell, 2012), human disturbance and distance to surface water (Odendaal-Holmes et al., 2014) play a role. Nonetheless, divergent views still prevail about the relative importance of factors underlying selective foraging of black rhinos and ultimately how they affect habitat choice (Muya & Oguge, 2000). We adopted a multiscalar approach to investigate black rhino forage selection and distribution across multiple scales (Figure 1). At the largest scale, we analyzed black rhino populations' spatial densities in relation to ecosystem productivity (Normalized Difference Vegetation INDEX [NDVI]). The species may be selecting highly productive areas in order to maximize biomass intake in minimal time (Bergman, Fryxell, Gates, & Fortin, 2001). Next, we examined black rhinos' distribution in relation to plant community selection. Specific plant communities may be selected over others for maximizing preferred food intake (Anderson et al., 2018). Then, we analyzed plant species selectivity by establishing seasonal diets through observational transects. Finally, at a molecular level we analyzed plant species' chemical profiles through metabolomics, because black rhinos, like other megaherbivores, may select specific plant species based on their primary and secondary chemical composition (Anderson et al., 2018;Muya & Oguge, 2000;Ndondo, Wilhelmi, Mabinya, & Brand, 2004). This dataset presents comprehensive information on black rhinos' distribution in relation to foraging across the hierarchical organization of the savanna ecosystem.

| Study site and vegetation communities
The study was conducted in Ithala Game Reserve (IGR), situated in Northern KwaZulu-Natal (27 30 0 S, 31 25 0 E) in South Africa with an area of 296 km 2 . Long-term annual rainfall averages 748 mm and the majority of this falls during the wet season (October-April), particularly, during the summer months (November-April) ( Figure S1). The reserve lies within the Savanna /Grassland biome and hosts 26 different types of vegetation communities (Van Rooyen & Van Rooyen, 2008) (Figure S2), mainly consisting of grasslands, thickets, rocky/ open/dense bushveld, woodlands, forests, riparian vegetation, cliffs, scarps, and open disturbed patches. Each individual rhino has been marked with a unique set of ear notches for identification and monitoring purposes, which allows the spatial position to be recorded every time an individual is sighted. We used black rhino monitoring data from 2012 to 2017 for this study. This sensitive data may not be published for the protection of the rhinos, however, the entire population on the reserve was monitored. This resulted in a total of 2,371 observation data points, with a median value of 47 observations per rhino for the 5 years of sightings. Black rhino location data was collected from dusk to dawn by all field rangers stationed across all quadrats of the reserve, thus all parts of the reserve were provided with equal monitoring effort. In particular, since IGR management's goal is to observe each rhino at least once a month, the monitoring effort across all individual rhinos and season should be uniform.

| Normalized difference vegetation INDEX (NDVI) calculation
As a general indicator of ecosystem productivity, or vegetation greenness, we used the Normalized Difference Vegetation INDEX (NDVI) calculated from Landsat TM, assuming that higher NDVI values indicate more productive areas (Berry, Mackey, & Brown, 2007). The NDVI was calculated using red (RED) and near-infrared (NIR) spectral bands of Landsat-7 with 30-m resolution: NDVI = (NIR − RED)/(NIR + RED). We employed Landsat data collected across the 5 years of rhino monitoring (2012-2017) to estimate a monthly average NDVI. Landsat data layers including more than 20% of cloud cover were excluded. We next used this dataset to characterize two distinct seasons; dry (May to October) and wet (November to April) ( Figure S3).

| Black rhino diet selectivity
To establish plant species selectivity, plants browsed and those avoided by black rhino individuals were surveyed on feeding paths (Shrader, Bell, Bertolli, & Ward, 2012). Feeding black rhinos were located by tracking them or through opportunistic sightings (Table S1). To limit autocorrelation, feeding paths from the same individuals were sampled at a minimum of 24 hours apart. Transects were approximately 50 m long and 2 m wide. Tracks determined the start and the F I G U R E 1 Multiscalar approach used in this study to explore factors driving black rhino habitat choice and forage selectivity. Each layer represents a scale of analysis applied to the study of the Savanna ecosystem direction of transects, predominately backtracking the animal's feeding path. A waypoint was recorded for each transect with the use of a handheld Garmin e-trex 20 GPS. All woody plants, shrubs and trees with a maximum canopy height of 2 m were recorded on the feeding path. Grasses were excluded as black rhinos are predominantly browsers (IUCN, 2011). Forbs were also excluded as we could not ascertain whether they were browsed by the rhinos when missing from the vegetation plot (Kotze & Zacharias, 1993). On the other hand, freshly browsed woody species are recognizable by the lighter color of the exposed wood, the wetness of remaining branches, and the characteristic way in which black rhinos browse. Black rhinos bite off large twigs in a pruning-shear manner attributable to the morphology of their hook shaped lips, leaving branches pruned at a clear 45 angle (Kotze & Zacharias, 1993;Shrader et al., 2012). Along each transect, all woody plant individuals were then scored as browsed or nonbrowsed for estimating plant species selectivity. Plant species selectivity was quantified by the application of Ivlev's electivity index E = (r i −p i )/(r i + p i ) (Strauss, 1979), where E is the measure of electivity per species, r i is the sum of browsed individuals of the same species and p i the relative abundance of the same species. The relative abundance of each species (p i ) was calculated by dividing the number of times it was encountered by the sum of encounters of all species. Selectivity is given and ranked for each species as the index has a possible range of −1 to +1, with negative values indicating avoidance of the plant species, zero indicating random selection from the environment and positive values indicating active selection (Strauss, 1979).

| Metabolomic analyses
To address chemically-based mechanisms of black rhino foraging decisions, we performed untargeted metabolomic analyses of the three most frequently browsed species (Lantana camara, Dichrostachys cinerea, and Acacia karroo, see Section 3), and three of the most frequently avoided species (Acacia gerrardii, Euclea crispa, and Lippia javanica). At least 10 leaves from five different individuals were collected for each species. All samples were dried and stored in brown paper bags before performing untargeted metabolomics analyses using UHPLC-QToF-MS as described in Data S1 Methods' section. On the same samples we also measured total carbon and nitrogen using a standard elemental analyzer (Flash 2000, Thermo Scientific, Waltham, Massachusetts, United States) to calculate the carbon to nitrogen ratio (C/N). First, we predicted that the metabolomics profile would differ between browsed and avoided species. Second, because black rhinos are deficient in nitrogen (as are most animals), we predicted that the browsed plants should contain higher levels of nitrogen content (thus a lower C/N ratio) than avoided plants.

| Statistical analyses
First, for each georeferenced rhino sighting, we estimated a specific point density based on the number of individuals observed in a 500 m radius (Wand, 1994) using the package pointdensityP (Evangelista & Beskow, 2018) in R-3.5.3 (R Development Core Team, 2018). Second, we calculated an index of selectivity for each plant community as the average of the species selectivity weighted by F I G U R E 2 Glmer output displaying; a, the relationship between black rhino density and NVDI; and b, the relationship between black rhino density and average vegetation type selectivity across seasons (yellow = dry and green = wet seasons) their relative abundance. The relative abundance of plant species in each plant community was estimated using transect data (number of stems divided by the cumulative length of transects per plant community). By doing so, we were able to extrapolate precise dietary information obtained during tracking at the local scale to the same plant communities' types. To estimate the influence of NDVI, selectivity and season on rhino density we fitted a Generalized Linear Mixed-Effects Models (glmer function in the package lme4 [Bates, Mächler, Bolker, & Walker, 2015]) with Poisson distribution as: density = NDVI*season + community selectivity*season + random/rhino identity. We visualized the output of the model using prediction of parameters interaction plots and river plots, which are particularly useful to represent interactions across large datasets (Hao, Dayal, Keim, Sharma, & Mehta, 2009), and illustrates the inner works of the model previously discussed. The column (d) shows how each vegetation type is constructed according to species selectivity and abundance. Each color represents the average species' selectivity for each species in a given community; blue = more attractive, and red = more repulsive. In the legend, "Aggregated" refers to locations that are more frequently occupied by black rhinos relative to the other locations, which are designed as "Dispersed." The metabolome of the six species sampled was analyzed using a bidirectional orthogonal partial least square discriminant analysis (OPLS-DA) with Pareto scaling (Thévenot, Roux, Xu, Ezan, & Junot, 2015) on the two foraging classes (browsed and avoided). OPLS has a similar predictive capacity compared to PLS by improving discrimination of the predictive components (Pinto, Trygg, & Gottfries, 2012). This model results in a two-dimension projection of plant chemical similarity, which allows visualizing rhinos' feeding choice in relation to plant metabolites identity. We used mixed-effect models to test whether rhinos' feeding choice (fixed factor: browsed versus nonbrowsed) was related to C/N content in plants. Plant species were included in the model as random factor (lmer function in the package lmer4 in R-3.5.3 [Bates et al., 2015]).

| Rhino density across NDVI and community selectivity
The spatial distribution analysis showed that rhino density is strongly variable and structured across the landscape (Figure 2, Figure 3a). In other words, across the reserve we observed locations that are more frequently occupied ("Aggregated" section in Figure 3a) relative to the others ("Dispersed" section of Figure 3a). The output of the mixed-effect model showed that black rhino density was negatively correlated with ecosystem productivity (NDVI) and positively correlated with plant community selectivity (Figure 2, Table 1). This was also illustrated by the river plot's ( Figure 3) flow that suggests that areas with higher density of black rhinos are correlated to low NDVI (Figure 3b), and plant communities with high global selectivity (Figure 3c). The opposite is true when observing the flow starting from the "Dispersed" rhino abundances.
We also found that community selectivity is driven by the balance between highly preferentially browsed (colder colors boxes), nonspecific (warmer color boxes) and avoided (red boxes) plant species (Figure 3d). Finally, we observed an effect of season on black rhino density,  where in the wet season black rhinos were present in higher numbers in low productivity areas and in plant communities with high selectivity values ( Figure 2, Table 1) than in the dry season. The wet season corresponds to a higher global productivity ( Figure S3).

| Species-specific forage selectivity
Over the course of two distinct seasons, black rhinos' diet proportions varied significantly ( Figure S4; Chisquare = 158.95, df = 70, p < .001), but overall, Ivlev's values for selectivity indicated that black rhinos preferentially browsed on certain plant species ( Figure S5). These The metabolomics results of tested samples (browsed and nonbrowsed) were plotted by chemical profile dissimilarity with an OPLS-DA. The latter highlighted a chemical discrimination between browsed and nonbrowsed species showing two clear clusters (Figure 4a), indicating that black rhinos select plants to forage based on specific metabolic profiles. On the other hand, while browsed plants have a 70% lower C/N compared to avoided plants (Figure 4b), we found no significant difference between browsed and avoided species across the six species tested (mixed-effect model for testing browsing preference effect; F 1,24 = 1.41, p = .37).

| Rhino spatial distribution
Through spatial analyses and the river plot, we were able to show that low densities of black rhinos were observed in vast areas that correspond to low selectivity plant communities. These findings concord with the fact that black rhinos are solitary sedentary animals and live alone within home ranges (Burt, 1943;Mitchell & Powell, 2012). Accordingly, to meet the nutritional requirements in low selectivity plant communities, they should maintain large foraging areas (le Roex, Dreyer, Viljoen, Hofmeyr, & Ferreira, 2019; Reid, Slotow, Howison, & Balfour, 2007). Additionally, we found that high-densities of black rhinos were associated with smaller areas of land and also with plant communities with high selectivity values. This may be because in these areas, nutritional requirements are met faster and therefore they can tolerate sharing their resources with other black rhinos. This effect is particularly marked in the wet season, indicating that during the dry season the decrease of plant biomass could be a limiting factor for rhinos to group more densely.

| Effect of productivity and selectivity on black rhino distribution
We found that black rhino density was negatively correlated with NDVI. This goes along with the observations that black rhinos favor open woodland and shrubland for optimal browsing as opposed to closed canopy woodlands (Gadiye & Koskei, 2016;Kotze & Zacharias, 1993). Open woodlands have lower NDVI values due to a lower presence of leafy plants and is therefore consistent to our findings of black rhinos' higher distribution in low NDVI areas. The species has also been reported to utilize dense bush for bedding sites (Anderson et al., 2018), which justifies moderate utilization of other areas. That said, if grassland NDVI values were to be excluded from NDVI plot averages, it is questionable whether black rhinos' choice of plant biomass would still qualify as low. Confirmation would require more intense tracking regimes across the year, and more accurate measurements of plant productivity at each site.
Black rhinos not only choose plant communities with low productivity but also high selectivity values. This finding is likely explained by the fact that they are selective feeders (Ganqa, Scogings, & Raats, 2005;Muya & Oguge, 2000), and the general fact that spatial variation in the availability of different food plants drives herbivore distribution (Kartzinel et al., 2015). Furthermore, a diverse and selective diet may be sustained by two constraints that black rhinos face. Firstly, being bulk feeders tolerant of low quality foods (Owen-Smith, 1992), a diverse diet is necessary for required nutrient intake (Muya & Oguge, 2000). Secondly, because Perissodactyls (odd-toed ungulates which include all rhinos) are unable to benefit from bacterial degradation of toxins (Freeland & Janzen, 1974;Muya & Oguge, 2000), a more diverse diet reduces the effect of ingesting high doses of toxic phytotoxins (Rhoades & Cates, 1976). We also found that the importance of plant community selectivity decreased during the dry season, suggesting that black rhinos must be less selective to cover nutritional requirements during this harsher season. All together this suggests that the ratio between preferred and avoided species determines the strength of community selectivity.

| Metabolomics and choice of forage
Previous studies (Ndondo et al., 2004;Van Lieverloo et al., 2009) indicated that black rhinos do not select forage to maximize nutrient intake and/or digestibility. Nutrient intake, measured by C/N ratio, has classically been related to feeding choice in invertebrate herbivores (White, 1984), and to a lesser extent for megaherbivores (Hopcraft, Anderson, Pérez-Vila, Mayemba, & Olff, 2012). Therefore, secondary metabolites are proposed as possible mechanisms of selection for black rhinos (Anderson et al., 2018). Correspondingly, our preliminary results suggest that black rhinos select their diet through chemistry rather than nutrient uptake. Secondary metabolites, such as phenol-based compounds, seem to play a bigger role than C/N ratio. Black rhino, in largely abandoning good vision capacities, have evolved highly sophisticated olfactory senses. Therefore, host plant choice should be largely mediated by odorous cues, but this hypothesis needs further investigation as well as to determine which compounds drive selection or avoidance.

| Study limitations and further steps
A problem encountered by this study, is the limited information on relative abundance of species in vegetation communities available for analysis. To deal with this, abundances on transects were projected across hosting vegetation types. But, microhabitats within vegetation types can be highly variable (Kartzinel et al., 2015) and transect locations were selected according to the presence of black rhinos. This might have led to a misrepresentation of microhabitat selection within a vegetation community. To counter this limitation, a meticulous and exhaustive vegetation survey covering the entire reserve is necessary. Some studies (Hall-Martin, Erasmus, & Botha, 1982;Oloo, Brett, & Young, 1994) imply that forbs could include an important fraction of black rhino diet. However, because of the unrecognizable way forbs are browsed due to their herbaceous nature (Kotze & Zacharias, 1993), they were excluded from surveyed vegetation transects. This may have impacted our results on diet proportions during the wet season, as there was a substantial presence of forbs on some of the transects. To resolve this issue and have supplementary accuracy in dietary information, dung samples could be analyzed in conjunction with vegetation surveys. Partially digested plant species fragments can be identified through microhistological examination (Van Lieverloo et al., 2009) or DNA analysis (Pompanon et al., 2012). A recent study (Anderson et al., 2018) showed that DNA analysis and forage plots results correlated well.
As a further matter, black rhino location data is based on sightings collected by rangers on patrol and its use is limited by two constraints. Firstly, because rangers are unable to identify individuals at night, data was usually collected during the day. Secondly, rhinos are less visible when the surrounding vegetation is thick (Walpole, 2002) and in spite of providing equal effort to all areas of the reserve, this could lead to occasionally overlooking rhinos. Black rhinos are diel (Joubert, 1971) and therefore collecting data during both day and night would not only enable finer scale modelling, but also show that feeding preferences, as their behavior (le Roex et al., 2019) might substantially differ across day and night. The use of GPS telemetry would greatly improve the temporality and accuracy of the dataset.

| Conservation implications
One of the key results of this study is that black rhino forage selection is driven by specific plant community selectivity and not productivity; and secondary compounds seem to play an important role in plant species selectivity. These findings must be incorporated into black rhino population management strategies and warrant further investigation. By considering precise black rhino habitat use across a hierarchical organization, carrying capacity calculations could alter. The identification of specific plants that black rhinos use preferentially, should assist in the development of more accurate reserve management programs that take into account plant community selectivity and productivity measures.
AUTHORS' CONTRIBUTIONS S. R., R. v. d. W., and V. D. initiated the project. V. D. collected the data. G. G. performed chemical analyses. E. M. and V. D. analyzed the data. S. R. and V. D. wrote the first version of the manuscript and all authors contributed to the final writing.

ETHICS STATEMENT
All data were collected following Ezemvelo KZN Wildlife research ethics. No ethical approval was required for the study.