Rolling pits of Hartmann’s mountain zebra (Zebra equus hartmannae) increase vegetation diversity and landscape heterogeneity in the Pre‐Namib

Abstract Microsites created by soil‐disturbing animals are important landscape elements in arid environments. In the Pre‐Namib, dust‐bathing behavior of the near‐endemic Hartmann's mountain zebra creates unique rolling pits that persist in the landscape. However, the ecohydrological characteristics and the effects of those microsites on the vegetation and on organisms of higher trophic levels are still unknown. In our study, we characterized the soil grain size composition and infiltration properties of rolling pits and reference sites and recorded vegetation and arthropod assemblages during the rainy season of five consecutive years with different amounts of seasonal rainfall. We further used the excess green vegetation index derived from drone imagery to demonstrate the different green up and wilting of pits and references after a rainfall event. In contrast to the surrounding grassland, rolling pits had finer soil with higher nutrient content, collected runoff, showed a higher infiltration, and kept soil moisture longer. Vegetation in the rolling pits was denser, dominated by annual forbs and remained green for longer periods. The denser vegetation resulted in a slightly higher activity density of herbivorous arthropods, which in turn increased the activity density of omnivorous and predatory arthropods. In times of drought, the rolling pits could act as safe sites and refuges for forbs and arthropods. With their rolling pits, Hartmann's mountain zebras act as ecosystem engineers, contributing to the diversity of forb communities and heterogeneity of the landscape in the Pre‐Namib.


| INTRODUC TI ON
Biopedturbation, the disturbance of soils by animals, is an important and often essential functional component of many ecosystems worldwide (Coggan et al., 2018). It determines the spatiotemporal characteristics of soil patches and thereby contributes to the ecosystems' diversity and heterogeneity (Mallen-Cooper et al., 2019). Examples of the ecological importance of biopedturbation can be found in all climatic zones, but the focus is on arid regions (Coggan et al., 2018;Mallen-Cooper et al., 2019). The spectrum of animals for these disturbances spans almost all animal groups (Coggan et al., 2018;Whitford & Kay, 1999) and even includes logical heterogeneity of the landscape (Bestelmeyer et al., 2006;Tongway & Ludwig, 2005). Small-scale topographical features such as depressions and changes in soil properties form favorable microsites with higher nutrient and water availability (Bestelmeyer et al., 2006) that affect density and composition of the vegetation.
The depressions created by animals collect runoff after rain and have better infiltration properties, increased soil moisture, and a higher nutrient content due to the accumulation of dung and urine (Mallen-Cooper et al., 2019). The removal of existing vegetation further changes the competitive balance of plant communities (Mallen-Cooper et al., 2019;Romero et al., 2015). Hence, such microsites offer more favorable conditions for plant recruitment and establishment than the surrounding matrix (Dean & Milton, 1991a, 1991bMallen-Cooper et al., 2019).
The vegetation of such microsites is therefore often characterized by higher annual plant species richness, higher biomass production, and different plant community composition (Mallen-Cooper et al., 2019). The more favorable conditions provided by these depressions are particularly important for annual forbs (Mallen-Cooper et al., 2019;Whitford & Kay, 1999), which account for the largest share of species richness in arid grasslands, contribute to soil organic matter, and provide shelter and forage for arthropods, as well as larger animals and mammals (Siebert & Dreber, 2019).
This in turn increases the activity density and diversity of arthropods (Nickell et al., 2018;Ruttan et al., 2016) or mammals (Ewacha et al., 2016). Such ecological effects can persist over many years or even decades (Coggan et al., 2018;Hastings et al., 2007;McMillan et al., 2011). Consequently, many soil-disturbing animals act as ecosystem engineers and play an essential role in maintaining the diversity and productivity of the respective ecosystem (Mallen-Cooper et al., 2019). Furthermore, such ecosystem engineers are often of particular concern for nature conservation (Hastings et al., 2007;Mullan Crain & Bertness, 2006), as their loss can lead to a deterioration in ecosystem function (Eldridge & James, 2009). Coggan et al. (2018) list 121 species in their ecosystem engineering database, of which many occur in arid or dry regions. A species not yet included is the Hartmann's mountain zebra (Equus zebra hartmannae), a species that is endemic in the semidesert of Namibia's great escarpment, the Pre-Namib, and a small part of Southern Angola and South Africa (Gosling et al., 2019). With their unique dust-bathing behavior, the Hartmann's mountain zebras create characteristic small depressions, which even when abandoned remain in the landscape for many years (Skinner & Chimimba, 2005).
The exact effects of such rolling pits on the vegetation of arid savanna ecosystems and within the trophic cascade are largely unknown and have not been reported yet. With our study, we aim to elucidate the role of Hartmann's mountain zebras as ecosystem engineers in arid savanna ecosystems. Therefore, we investigated the effect of the rolling pits of Hartmann's mountain zebras over a period of five consecutive vegetation growing seasons with different amounts of rainfall on soil properties and density and composition of the vegetation. Further, we tested whether and how these soil and vegetation differences also affect ground-dwelling arthropods.

| Study site and study species
Our study was conducted between 2014 and 2018 on the farm Rooiklip, Khomas,Namibia (33 K 612,448 7,411,832), situated at 1,000 m.a.s.l. in the center of the escarpment halfway between Windhoek and Walvis Bay ( Figure A1). Namibia's great escarpment is an approximately 100-km wide strip of rugged semidesert reaching from the border with South Africa in the South to Angola in the North. The climate is hot-arid, mean annual precipitation ranges from ~50 mm in the West to ~200 mm at the upper edge in the East (Mendelson et al., 2009), but rainfall is highly variable and often erratic. At our study site, the mean annual precipitation is 120 mm, but ranged from 60 mm to 180 mm during the study period.
Rainfall was measured on a daily basis with an on-site standard rain gauge. The main vegetation growth period occurs from February to March, where two-thirds of the annual precipitation occurs (Wagner et al., 2016). The vegetation is arid grassland dominated by perennial grasses of the genus Stipagrostis that is sparsely scattered with trees and occasional dwarf shrubs. Depending on the amount of rainfall, various annual forbs and annual grasses join these perennials. Total vegetation cover seldom reaches more than 15% (Wagner et al., 2016). Soils are leptosols and consist of a shallow layer of differently sized weathered particles of the underlying schist bedrock.
Rooiklip and its neighboring farms are largely unfenced and allow wildlife to roam freely between the farms and from the edge of the escarpment to the Namib Desert.
The escarpment and the Pre-Namib are home to the near-endemic Hartmann's mountain zebra ( Figure A1; Environmental Information Service, 2020; Gosling et al., 2019). The habitat of Hartmann's mountain zebra is mountainous areas and flats. Hartmann's mountain zebras further exhibit a marked migration pattern, following the rainfall and forage availability and moving between dry season and rainy season home ranges . Hartmann's mountain zebras usually live as family groups that consist of a stallion and his harem of mares and their foals. However, stallion groups or solitary stallions are frequent (Skinner & Chimimba, 2005), and occasionally, several family groups can join into herds with up to 30 individuals (Skinner & Chimimba, 2005).
Similar to many other equines, the Hartmann's mountain zebras show a regular dust-bathing behavior (Joubert, 1972a), that in contrast to many other equines is carried out quite frequently several times during the day: The animal first moves into an incomplete lateral position and then turns alternately around its longitudinal axis to the left and right and then gives itself a push and stretches itself (McGreevy, 2012;Panzera et al., 2020). The exact reason for this behavior is unknown, but it is assumed that it is mainly for selfgrooming and removal of ectoparasites but partially also for resting (Joubert, 1972a). Typically, a family of Hartmann's mountain zebras prefers a very specific area for this procedure and often the individual spots are used by the group members, so that soon pronounced depressions in the soil are created that have on roughly 2.5 m diameter and can become up to 30 cm deep (Joubert, 1972b). When the zebras are migrating, the pits of the previous home range become abandoned but are sometimes reactivated during the next season, where they are used by the same or other families.

| Microsite characteristics
In 2014, we mapped all detectable zebra pits within a 1 × 1 km area unused for livestock-keeping since 2001. We searched the pits for signs of current use (tracks, fresh dung, lack of vegetation), noted the presence of zebra feces, and classified the pits as either active or abandoned. From all pits mapped in 2014, we randomly selected 16 abandoned pits, marked them, and paired them with a 2 m × 2 m reference site 5 m south of the center of the respective pit ( Figure A4).
We chose this particular size to match the average size of the rolling pits. The distance of 5 m was chosen to ensure similar topography and ground conditions. We measured diameter and depth of the depression and determined the depth of the loose soil layer of the pit and reference sites using a simple penetrometer (Sutherland, 2006). We further characterized soil grain size composition of pits and reference sites by taking soil samples of the upper ten-cm soil layer of a 40 × 40-cm-sized area. The samples were sieved into the different grain size classes: cobble and boulders, very coarse gravel, coarse gravel, fine-tomedium gravel, and sand and smaller particles, according to ISO 14688 (International Organization for Standardization, 2017), and the respective percentage by weight was calculated. Water-soluble soil nutrients were determined by eluting an aliquot of 10 g the soil fraction sand and finer particles (oven-dried at 105°C for 24 hr) with 100 ml of ultrapure water for 3 hr at room temperature. The eluate was left to sediment for 20 min, the supernatant was centrifuged for 15 min at 3,500 g, and 5 ml of the supernatant was used to analyze nitrogen and phosphate contents using standard ion chromatography (Michalski et al., 2019). A simulated precipitation event equivalent to 20 mm rainfall (applied within 1 hr) was used to measure the infiltration depth and determine the time course of soil moisture within the upper 10 cm of soil. Soil infiltration depth was determined 30 min after the precipitation event by digging and measuring the depth of the visually identified seepage front. The time course of soil moisture was determined using a TDR probe (Delta-T Devices, ML3; Rajkai & Ryden, 1992) in two separate runs over 6 days each. The probe was initially calibrated once with an aliquot of soils taken from the pits and reference sites, respectively (Delta-T Devices Ltd, 2013).

| Vegetation assessment
From 2014 until 2018, we consecutively mapped the vegetation of pits and reference sites in two rounds at an interval of about 4 weeks during the respective vegetation period (Table A1). We determined the plant species present and their respective abundance and visually estimated the percentage of cover. Cover data were aggregated to the percentage of total plant cover, cover of annual forbs, and cover of annual grasses.
To show the development of vegetation over time and the differences in greening and wilting of pits and reference sites, we used high-resolution RGB images acquired with a commercial DJI Phantom 4pro UAV. During the rainy season 2017, between 27 February and 6 April 2017, in total, 15 flyovers, 2 to 3 days between each other, were carried out over a transect covering the study area. The last relevant rain event occurred 14 days before the start of the flyovers and was 8 mm. During the flyovers, four distinct rainfall events occurred with 5 mm on 28 February, 22 mm on 2 March, 9 mm on 8 March, and 4 mm on 13 March. Flight altitude was 40 m.a.g.l., and image overlap was 75% to allow proper image alignment. All images of each flyover were processed and stitched as described in Woellner and Wagner (2019), resulting in a georeferenced (UTM, 33S) orthomosaic with a resolution of ~2 px/cm and a digital elevation model (DEM) of the study area.

| Arthropod sampling
Coincident in time with each vegetation mapping, the occurrence, and abundance of ground-dwelling arthropods were determined using pitfall traps. As the numbers of trapped arthropods depend not only on their abundance but also on their activity levels (e.g., Topping & Sunderland, 1992), the obtained numbers of arthropods are commonly referred to as "activity density" rather than abundance (Thiele, 1977). The resulting data were used to test for differences among arthropod assemblages between pits and references and to study bottom-up effects along the trophic cascade.
Pitfall traps had a volume of 650 ml and a diameter of 9 cm and were filled with 250 ml of a 50% ethylene glycol-water mixture and a drop of detergent to reduce surface tension. Traps were installed on a level with the soil surface in the center of each pit or reference site (Lange et al., 2011) and kept open for 7 days after the vegetation was mapped. After each sampling round, arthropods were transferred into 70% ethanol and identified up to order and classified according to feeding guilds into herbivores, omnivores, and predators (Table A2; Picker et al., 2004;Scholtz & Holm, 2008) following the rapid assessment of biodiversity proposed by Oliver and Beattie (1996) used in comparable studies, for example, by Nickell et al. 2018 on Bison wallows. Though this somewhat limits deductions in the community structure, this was necessary as the majority of invertebrates have not been described (Samways & Samways, 1994).
However, the method has been successfully applied to demonstrate habitat-related changes in arthropod assemblages (see Blaum et al., 2009;Fabricius et al., 2003). Beetles (Coleoptera) and true bugs (Hemiptera) that included both species with herbivorous and predatory feeding preference were assigned to omnivores. Larvae that were not further determinable were excluded from analysis.

| Statistics and data analysis
Statistical analysis was performed with R, version 3.6.3 (R Core Team, 2020). Comparisons of soil grain sizes and nutrients between pits and reference sites were done by permutational paired t tests (perm.t.test, library RVAideMemoire, Hervé, 2020). Time course of soil moisture was modeled as linear mixed effect model (lme; library nlme version 3.1-120; Pinheiro et al., 2020). Soil moisture was log(x+1) transformed to obtain normality of variances. Day (since simulated rain event) and Type (pit/reference site) were used as explanatory variable.
Plant diversity was calculated as "effective number of species" (ENS), the equally common species that result in the respective Shannon-Wiener index and reflect the true diversity of the communities considered (Jost, 2006) using the function diversity from the vegan library (Oksanen et al., 2019). Indicator species, that is, plant species predominantly associated with either pits or reference sites, were identified using the multipatt function with IndVal.g as association function (library indicspecies; De Caceres & Legendre, 2009).
The function provides estimates for the probability of a species occurring in the respective site group (pO) and a probability with which a site belongs to the respective group when the species is found (pT). It further calculates for each species an indicator value (IndVal) based on the product of its relative frequency and its relative average abundance in either pits or reference site.
Analysis of the drone data was done using the raster package (library raster, version 3.1-5; Hijmans, 2020). Based on the orthomosaic, the "greenness" of each pixel was determined by using the calc function to calculate the excess green minus excess red vegetation index (ExGR; Meyer & Neto, 2008) that is suitable to differentiate Differences between greenness of pits and reference sites were tested for each day using permutational t tests with holm adjustment for paired samples; differences between greenness of pits or ref- erence sites to baseline were tested using standard permutational t tests (perm.t.test, library RVAideMemoire, Hervé, 2020) with holm adjustment.
To characterize the vegetation composition with regard to the prevalence of a certain functional group (annual grass, annual forb), we calculated the annual grass-to-forb ratio for each pit and reference site as the difference between the cover of annual grasses (annual grass cover/total cover) and proportional forbs cover (annual forb cover/total cover). This resulted in an index ranging between −1 for forbs only and +1 for grasses only. The relation between this grass-to-forb ratio and rainfall received within the preceding 30 days was modeled as linear mixed effect models (lme; library nlme version 3.1-120; Pinheiro et al., 2020) using maximum likelihood with amount of rainfall (Rain30d), Type, and the two-way interaction as explanatory variables.
For testing bottom-up effects on arthropods, the activity density of the respective feeding group was modeled as linear mixed effect model (lme; library nlme version 3.1-120; Pinheiro et al., 2020) using maximum likelihood. Models for herbivores included total vegetation cover and Type (pit/reference site) as explanatory variables, models for omnivores additionally contained the activity density of herbivores, and models for predators additionally contained the activity density of potential prey (herbivores and omnivores), as well as two-way interactions of all explanatory variables. To obtain normality of errors, the response variables were log(x+1) transformed.
For all lmes, we used maximum likelihood and random intercept with pit/reference id as random factor to ensure independence of errors with respect to spatial and temporal autocorrelations (Pinheiro & Bates, 2000). Model simplifications were done in a backward stepwise model selection procedure by stepAIC (library MASS; Venables & Ripley, 2002) until the respective minimal adequate model was reached.

| RE SULTS
Between 2014 and 2018, in total, 656 pits were mapped in the study area, all restricted to elevated, but level and flat grounds, free of rocks, with an average inclination of <5° ( Figure A5). From these 656 pits, 341 (~52%) were created in new places, where no pit has been mapped before, the rest was either re-used or abandoned during the study period. Fresh or dried dung was found in 89% of the mapped pits. After heavy rainfall events, many rolling pits collected runoff and rainfall and it took up to three days for the water to evaporate ( Figure A3b).
Pit diameter ranged from 105 cm to 350 cm with a mean of 229 ± 55 cm (mean ± SD), resulting in an average area of 4.1 ± 0.3 m² per pit. The average depth of the pits (distance from the surface of the matrix to the soil layer of the pit) was 9.5 ± 3.6 cm. The substrate of pits contained considerably more fine-grained material than the reference sites. In pits, sand and finer grained particles (84%) were the highest share, followed by gravel (11%), whereas reference sites contained 60% gravel and only 36% sand and finer material (Table A3). Water-soluble soil nutrients were low (Table A4), with nitrogen and phosphate contents being significantly higher in pits (water-soluble N: 4.8 mg/kg; PO 4 : 1.3 mg/kg) than on reference sites (water-soluble N: 1.2 mg/kg; PO 4 : 0.5 mg/kg). The depth of the loose soil layer of pits was 29.3 ± 8.1 cm and thereby significantly higher than the average loose soil depth 5.5 ± 2.1 cm of the reference sites (p < .002). After our simulated rain event, the soil infiltration depth in pits reached 9.9 ± 1.9 cm and only 4.7 ± 1.7 cm on reference sites, which corresponded with the depth of loose soil. During the first 24 hr after the event, the average soil moisture in pits reached 6.31 ± 0.24% and remained above the permanent wilting point for coarse sand (~2%) for over four days. In contrast, soil moisture never exceeded the permanent wilting point on reference sites (Figure 1).
The temporal course of greening and drying out of the vegetation following natural rain events showed clear differences between pits and reference sites. The ExGR vegetation index and hence greenness of the zebra pits increased with each rainfall event until day 14 and then decreased steadily over the next 14 days until it reached the baseline at day 28 when the vegetation was wilted.
From day six on, it was always above the greenness of the reference sites. In contrast, the greenness of the reference sites increased only after the first rain and then remained constant, apart from a slight decrease on days six and eight, until day 18. After that, it decreases F I G U R E 1 Trend of soil moisture in pits and reference sites after a simulated rainfall event with 20 mm applied within 1 h. Dashed horizontal line indicates permanent wilting point of sandy "soils"  Table A5).
Total vegetation cover did not exceed 23% in pits and 13% on reference sites and correlated with rainfall received within the preceding 30 days. Both total vegetation cover and cover of annual forbs were higher (total cover on average 1.4-fold, cover of annuals 1.8-fold) in pits compared to reference sites throughout all years (Table A6). The proportion of forbs and the total cover of annual plants was significantly higher in pits and increased significantly with the amount of rainfall received during the preceding 30 days ( Figure 3; Table A6).
In total, we found 26 annual plant species, from which six were annual grasses and 20 annual forbs, but all native to the area. General plant diversity of pits and reference sites was similar with 4.8-8.4 species in pits and 4.5-7.1 species in reference sites (Table A6) but varied throughout the years as a response to different amounts of rainfall. Pits were mainly characterized by the presence of annual forbs (seven indicator species), whereas annual grasses mainly characterized the reference sites (three indicator species; Table 1).
Herbivorous arthropods showed a significant positive correlation with total vegetation cover with slightly higher activity density in pits compared to reference sites. Omnivores were positively related to the activity density of herbivores, but without differences between pits and reference sites, whereas predators were showing a significant positive correlation with the activity density of prey (omnivores and herbivores) but no significant difference between pits and reference sites (Table 2).

| D ISCUSS I ON
The characteristic rolling pits of the Hartmann's mountain zebra in the Pre-Namib, created by their dust-bathing behavior, produce favorable microsites that clearly differ from their surroundings.
Their special ecohydrological properties promote a denser and different vegetation composition, characterized by annual forbs, and, to a lesser extent, increase the activity density of herbivorous ar- The characteristic of the rolling pits generally corresponds to those known from depressions and resting places of larger antelopes, equines, or kangaroos. Similar to the bedding sites of the Nubian ibex (Gutterman, 1997) or the wallows of American bison (McMillan et al., 2011), the rolling pits of Hartmann's mountain zebra collect runoff after rainfall. The loose soil filling of the pits allows a better infiltration, and once infiltrated, this moist soil acts as a storage for the runoff water that keeps the soil moist for a few days longer than the surrounding soil. In addition, the loose soil of rolling pits has a higher water-soluble nitrogen and phosphate contents.
Nine of ten pits we mapped had considerable remains of dry zebra dung. Deposits of feces and higher urine concentrations in the soil are a typical feature known from other resting places of larger animals such as Gemsbok (Dean & Milton, 1991a, 1991b or Kangaroo (Eldridge & Rath, 2002). In addition, such depressions frequently accumulate litter (Mallen-Cooper et al., 2019), which is then together with the dung quickly converted into nutrients by the stronger mineralization due to a higher soil moisture (Eldridge & Rath, 2002;Veldhuis et al., 2018). As a result, these resting places and so the rolling pits of the Hartmann's mountain zebra form relatively moist and nutrient-rich patches in the otherwise nutrient-poor savanna (Augustine & McNaughton, 2006;McNaughton et al., 1988). The general importance of animal-created depressions in arid environments for annual forbs has been documented for a number of animal species (Coggan et al., 2018;Gutterman, 2001). The higher and longer lasting soil moisture of the rolling pits also plays a crucial role in the occurrence of annual forbs. Together with higher nutrient availability, it is a decisive factor for the establishment of forbs in arid environments (Siebert & Dreber, 2019), which need longer lasting favorable conditions for the successful establishment of their seedlings (Hillel & Kozlovsky, 2012). The lack of perennial grasses in the pits further ensures that freshly germinated ephemerals are subject to less competition (Dean & Milton, 1991a, 1991b. Many of the annual forbs found in the pits such as H. ciliatum, H. modesta, or  Back-transformed, exponentiated coefficients, and levels of significance are given. "n.t.": not tested. "Herbivores" include the orders Caelifera, Cicadina, and Stenorrhyncha; "omnivores" include Blattodea, Coleoptera, Ensifera, Heteroptera, and Psocoptera; "predators" include Acarina, Araneae, Chilopoda, Pseudoscorpiones, and Scorpiones.
*p < .05, **p < .01, ***p < .001. a "pit" versus "reference site"; reference level was "reference site." as an ecosystem engineer. In the future, the importance of the rolling pits for the local forb and arthropod communities might even increase, due to prolonged severe drought periods, which are expected in the region as a result of climate change (Archer et al., 2019;Maúre et al., 2018).

ACK N OWLED G M ENTS
We are grateful to the owners of Rooiklip, Hannelore Neuffer, and

CO N FLI C T O F I NTE R E S T
The corresponding author confirms on behalf of all authors that there have been no involvements that might raise the question of bias in the work reported or in the conclusions, implications, or opinions stated.  F I G U R E A 4 Aerial photograph with densely vegetated rolling pits. Details from an aerial photograph of the study site with zebra tracks (dotted lines) and six abandoned rolling pits (dashed circles) of mountain zebra 14 days after a rainfall event with ~30 mm of rainfall within 1 day. The dense annual vegetation clearly stands out against the less dense cover of the perennial Stipagrostis grasses of the matrix. The configuration of a research plot consisting of pit and reference site and the position of the pitfall traps (red dot) in the center are schematically illustrated F I G U R E A 5 Position of the rolling pits within a section of the study site. Rolling pits only occur on level ground and often elevated grounds that make up about 1/3rd of the study site. Ground with an inclination of less than 5° is highlighted with magenta, red dots mark the respective pits Note: Aliquot of 10 g oven-dried substrate (sand and finer) eluted with 100ml of ultrapure water for 3h at room temperature; eluate left to sediment for 20min, supernatant centrifuged for 15min at 3,500 g; 5m of supernatant analyzed using standard ion chromatography. Significance of difference between pits and reference determined with paired permutational t tests (perm.t.test, RVAideMemoire; n = 999); * p < .05, ** p < .01, *** p < .001; n.s.: not significant. Note: n.s.: not significant, ***p < .001; **p < .01, *p < .05. Average ExGR values and standard deviation of pits and references during the flyovers and the significance of difference between pits and references and the respective difference to baseline.