Effects of free‐ranging livestock on occurrence and interspecific interactions of a mammalian community

Abstract Mammalian communities inhabiting temperate grasslands are of conservation concern globally, especially in Central Asia, where livestock numbers have dramatically increased in recent decades, leading to overgrazing and land‐use change. Yet, how this pervasive presence of livestock herds affects the community of wild mammals remains largely unstudied. We used systematic camera trapping at 216 sites across remote, mountainous areas of the Mongolian Altai Mountains to assess the spatial and temporal patterns of occurrence and the interspecific relationships within a mammalian community that includes different categories of livestock. By adopting a recently proposed multispecies occupancy model that incorporates interspecific correlation in occupancy, we found several statistically strong correlations in occupancy among species pairs, with the majority involving livestock. The sign of such associations was markedly species‐dependent, with larger wild species of conservation concern, namely, snow leopard and Siberian ibex, avoiding livestock presence. As predicted, we found evidence of a positive correlation in occupancy between predators and their respective main prey. Contrary to our expectations, a number of intraguild species pairs also showed positive co‐occurrence, with no evidence of spatiotemporal niche partitioning. Overall, our study suggests that livestock encroaching into protected areas influences the whole local community of wild mammals. Though pastoralism has coexisted with wildlife for millennia in central Asian grasslands, our findings suggest that policies and practices to decrease the pressure of livestock husbandry on wildlife are needed, with special attention on large species, such as the snow leopard and its wild prey, which seem to be particularly sensitive to this pervasive livestock presence.


Introduction
Species are not distributed at random across space, but form communities whose structure and composition are moulded by biotic interactions (Davis et al., 2018;Wisz et al., 2013). With the dramatic increase of the human footprint on the planet (Venter et al., 2016), however, anthropogenic disturbance, particularly land use change, has become a prominent factor altering inter-specific interactions, and even community composition (e.g., Kiffner et al., 2014;Rovero et al., 2020;Suraci et al., 2021). Human activities such as hunting can for example affect mammalian communities through apex predator removal, causing the increase of large herbivore density and in turn a change in vegetation (Hebblewhite et al., 2005). Top-down perturbations can also lead to mesopredator release with higher pressure on small prey and other alterations of ecosystem functioning (Brook et al., 2012;Newsome et al., 2017). Another relevant anthropogenic activity is livestock grazing, as in some regions like central Asia it has increased steadily over the last decades, raising concern over the effects it may exert on communities of wild mammals (e.g., Berger et al., 2013). However, how grazing livestock affects the spatio-temporal occurrence and inter-specific interactions of whole communities of mammals remains largely unstudied.
Mammalian communities inhabiting temperate grasslands are indeed of particular conservation concern globally, as these constitute one of the most endangered and understudied biomes worldwide, threatened by land use change, overgrazing by livestock, and climate change (Hoekstra et al., 2005;Ritcher & Osborne, 2014;Nunez et al., 2020). Palearctic prairies are characterised by strong seasonal variations in precipitation and temperature, hence in grass greenness and palatability, driving herbivores' use of space (Mueller et al., 2008). On the other hand, pastures are influenced by herbivore grazing: for example, in central Asia pikas (Ochotona spp.) and Brandt's voles (Lasiopodomys brandtii) impact vegetation community structure and composition through burrowing (Bagchi et al., 2006;Sawamukai et al., 2012). The presence and density of herbivores in turn influence the distribution and abundance of carnivores (Ross et al., 2010a;Chetri et al., 2017). In these grasslands, livestock grazing is a major factor at play, especially in Central Asia where the number of domestic animals has dramatically increased in recent decades following the globalization in the cashmere wool market (Berger et al., 2013). Detrimental effects of livestock encroachment and competition on mammals have been documented, particularly on wild ungulates Ekernas et al., 2017) and large carnivores (Sharma et al., 2015).
However, the overall effect on the whole community of medium and large-sized mammals is poorly understood, and evidence shows that livestock may also alter the interactions among wild species (e.g. Rottstock et al. 2020;Feng et al. 2021).
In this study, we aimed to assess the spatial and temporal patterns of occurrence and the inter-specific relationships within a mammalian community in central Asia that includes livestock.
We sampled four areas in remote and rugged environments within the Mongolian Altai mountains through systematic camera trapping on a total of 216 sampling sites. We took advantage of the ability of camera traps to detect multiple species, including domestic ones, simultaneously and at the same spatial scale. We analysed data with a recently proposed multi-species occupancy model that accounts for imperfect detection and incorporates inter-specific correlation in occupancy (Tobler et al., 2019). This analytical approach allowed us to study the effects of both environmental and anthropogenic factors on occupancy at community-and species-level while quantifying the strength of potential interactions among species, as inferred from patterns of species co-occurrence (Richmond et al., 2010).
Our main research questions and predictions were: (1) Is there strong evidence for spatiotemporal interactions between wildlife and livestock? We predicted a generally negative cooccurrence of wild species with livestock. (2) If associations between livestock and wild species exist, are these related to species' body size? We predicted that larger species would display more negative associations with livestock. (3) Do predators and prey in the community co-occur in space and time? We expected a positive co-occurrence of predators and their respective prey. (4) Is there evidence of intra-guild competition through spatial and/or temporal avoidance? We predicted a negative co-occurrence in space and/or time of species of the same guild as a result of niche partitioning.

Study areas
This study was conducted in the Mongolian Altai mountains, in the western provinces of Bayan Olgii and Hovd, that are characterised by a mosaic of high plateaus and mountain ranges. The climate is semi-arid, with short summers and long, cold and dry winters. In the four study areas most inhabitants are nomadic herders. The pastoralists in the Sutai massif belong to the ethnic group of the Mongols, and profess Buddhism, while in the other areas they mainly belong to the ethnic group of the Kazakhs, professing Islam. In the sole Bayan Olgii province, where the first three surveyed areas occur, 2.2 million livestock were reported in 2018, 1.9 million of which were sheep and goats (National Statistics Office of Mongolia, 2019). Largesized livestock, represented by horses, camels, cattle and yaks, are mostly free ranging, while sheep and goats are guarded by herders and dogs during the day and held in corrals at night, except in summer (July-August), where usually no fence is used (Augugliaro et al., 2020). Mongolian National Parks are subdivided into three zones with different regulations: special zone, travel/tourism zone, and limited use zone. Traditional pastoralism is only allowed in the limited use zone, whereas it is not allowed at all in the Strictly Protected Areas.

Data collection: camera trapping and covariates
For each area we first created in GIS a regular grid with cells of 4 km 2 that aimed to cover at least the minimum extent of 500 km 2 recommended for snow leopard (Panthera uncia) population estimate studies (Jackson et al., 2005). Within each cell we located in the field one camera-trap site.
The number and locations of chosen sites were constrained by terrain morphology, snow depth and the number of camera-traps available for sampling (see Appendix S1: Table S1). Two of the PAs were too large to monitor their entirety, thus we selected a portion that we considered optimal habitat for the snow leopard according to previous research and local knowledge. Specifically, in Siilkhem-B we did not sample the southern end of the PA, that includes low elevation and less suitable habitat, and in Tavan Bogd we focused on the north-eastern part of the PA, since the southern part hosts mainly coniferous forests. In Khork Serkhe we sampled the whole PA and surrounding suitable areas, and in Sutai we sampled the entire mountain massif. Overall, the area sampled ranged in size from 513 to 1110 km 2 (Appendix S1: Table S1). As a result, we sampled with camera traps a total of 216 sites (48, 44, 63 and 61 for the four areas, respectively) that ran on average for 65 days each (SD = 20.7). The minimum distance between contiguous sites was 1.5 km.
We placed camera traps on small rock piles c. 50 cm above the ground at a distance of approximately 2-4 m from a target trail, setting no delay between consecutive triggers. Sites were usually located on passes, ridges, narrow passages or valley bottoms and were chosen to maximise the detection of snow leopards, which we considered to be the most elusive species. However, our data show that these sites were used by several medium to large sized mammals and by livestock.
We did not use any baits or lures. Camera-trap models were of three different manufacturers, characterized by high (Reconyx) and medium (Cuddeback and Bushnell) trigger speed. We annotated camera-trap images using the software Wild.ID (Fegraus & MacCarthy, 2016), and screened the data to derive for each area the checklist of species and their total number of detections. We set the temporal resolution of detections to 1 day, thus repeated detections of the same species on the same day were counted as a single detection event, as is common practice in occupancy modelling.
Using a Digital Elevation Model and a satellite map in the software Quantum GIS (QGIS, 2018) we derived three environmental covariates for each camera trap site: elevation, terrain slope and distance from the closest herders' settlement (consisting of gers, houses, villages or mining camps), whose geographic localization was acquired in the field with a GPS device.

Data analyses
Out of all species detected (Appendix S1: Table S2) we excluded those found at only one or two areas, which generally yielded <5 detections (namely Ursus arctos, Lynx lynx, Ovis ammon, Mustela erminea and Sciurus vulgaris), and focused on the following 12 species of medium to large sized mammals detected across at least three areas: one large herbivore (Capra sibirica), 3 large carnivores (Canis lupus, Panthera uncia, Gulo gulo), 4 medium carnivores (Martes foina, Mustela eversmanni, Otocolobus manul, Vulpes vulpes) and 4 medium to small herbivores (Lepus spp., Marmota spp., Ochotona spp., Spermophilus spp.). The last four were either represented by more than one species of the same genus occurring in different areas (i.e., Marmota sibirica was found in all areas except for Tavan Bogd NP where M. baibacina occurs, and Lepus timidus was found in the first two areas whereas L. tolai in the latter two) or occurring in the same area but undistinguishable from one another in camera-trap photographs. Hence, given the similar ecology of these congeneric species we merged them at the genus level. Moreover, we merged the detections of horses, yaks, camels and cattle into the category 'Large-sized livestock' and those of sheep and goats into 'Smallsized livestock'. The subdivision of domestic animals in these two classes is justified by previous research, which found that large and small-sized livestock are managed differently by herders and may therefore affect wild species' spatio-temporal patterns differently (Augugliaro et al., 2020).
We arranged the number of daily detections per species at camera-trap sites into a matrix of n = 216 sites by S = 14 species. We modelled these detections/non-detections using the joint species distribution model with species correlation and imperfect detection developed by Tobler et al. (2019). This framework is an extension of the multispecies occupancy model (Dorazio & Royle, 2005) that explicitly models the residual correlation in species' occupancy probability (c), corrected for imperfect detection. Our modelling approach is detailed in Appendix S2.
We tested the effect of study area, elevation, terrain slope and distance from herder settlements on occupancy, and of camera-trap sensitivity and distance from herder settlements on detection probability, both at community and species level. The chosen covariates were not correlated according to Spearman's correlation test (|rho| < 0.5). We then calculated the matrix of residual inter-specific correlations in occupancy with their corresponding BCI (Bayesian credible intervals). We implemented the model in a Bayesian framework using JAGS (Plummer, 2003) via R (version 3.6.2; R Development Core Team, 2019) with the package 'R2jags' (version 0.5.7-7; Su et al., 2015). We specified vague prior probabilities for occupancy and detection community hypermeans through normal distributions centred on zero, through uniform priors in the interval from 0 to 10 for the variance of species-specific random effects, and through uniform distributions in the interval -1 to 1 for the latent variables' regression coefficients (except for the mathematical constraints required for parameter identifiability, see Tobler et al., 2019). The regression coefficients of large-and small-sized livestock were not extracted from the community hyper-means, but were instead treated separately, since their inclusion could have biased the estimation of the response of the wildlife community to covariates. Community hyper-means are therefore to be interpreted as mean effects for the assemblage of wild species only. We generated three parallel chains of 200,000 iterations with a burn-in of 20,000 iterations and thinning by 10 to derive summaries of parameter posterior distribution. We checked for convergence of the Markov chains based on the Gelman-Rubin statistic (Gelman et al., 1992) and with visual examination of the chains.
To explore the potential relationship between species-specific body mass and the correlation in occupancy with livestock, we regressed the estimated correlation coefficient of each wild species with the two livestock classes against the log of the species' body mass, which we sourced in Smith et al. (2003) and Hayssen (2008).
We also analysed the temporal activity pattern of each species using a non-parametric Kernel Density Estimation function implemented with the package 'overlap' (Ridout & Linkie, 2009) in R. We first extracted independent events by aggregating detections separated by less than 30 minutes, and then created the activity distribution curve of each species. To evaluate potential patterns of temporal avoidance or coexistence between the species in the target community we computed the coefficient of overlap for each pair of species. The coefficient of overlap Δ ranges from 0 (no overlap) to 1 (complete overlap) and is obtained performing pairwise comparisons of diel activity patterns.

Results
Out of the 91 estimated pairwise residual correlation coefficients in occupancy (c), 24 had a 90% BCI that did not overlap zero (Figure 2). Of the correlations in occupancy between wildlife and livestock of both classes, 12 had a 90% BCI that did not overlap zero (see Figure 2  The correlation in occupancy between wild species and livestock appeared negatively related to species body mass, a pattern that was more marked when considering small-sized livestock, as illustrated in Figure 4. Hence, larger wild species tended to have a more negative correlation in occupancy with large-and especially small-sized livestock.

Discussion
By using camera-trapping data, activity pattern analysis and a multi-species occupancy model that incorporates inter-specific correlation in occupancy and imperfect detection, we studied the spatiotemporal patterns of occurrence and inter-specific relationships within a mammalian community in the Mongolian Altai.
Each of the 14 species and livestock categories considered showed strong statistical evidence for spatial co-occurrence (i.e., 90% BCI of c that did not overlap zero) with at least another species, with the only exception of the steppe polecat (Mustela eversmanni). Remarkably, of the 13 interactions that each species can display with any other, the largest number of statistically supported ones involved one of the livestock categories (8 and 6 for large-and small-sized livestock, respectively), while the highest number of interactions involving wild species was five, for both snow leopard and marmot. Considering such scores as a simple measure of the magnitude of spatial co-occurrence among species regardless of their signs, the results support our choice of including livestock in the analysed community, given their widespread presence in the region and the potential influence on wildlife (e.g., Berger et al., 2013;Soofi et al., 2018).
Large areas of Mongolia have been subject to semi-nomadic pastoralism for centuries, but the recent increase in livestock numbers has raised concerns about possible degradation effects (Wesche et al., 2010). While the numbers of cattle and camels remained relatively constant across the country, a steep increase in sheep and goat populations has been recorded since the early 1990s, with the latter likely to represent a greater impact to wildlife than large-sized livestock given their much larger herd size (Hilker et al., 2014). Indeed, we found that the relationship between residual correlation in occupancy of wild and domestic species and wildlife body mass was negative, with a more pronounced effect for the small-sized livestock. This pattern points to a higher sensitivity of larger species to livestock encroachment into protected areas, which is of particular concern for the snow leopard and the Siberian ibex, and more generally mirrors known patterns of size-dependent vulnerability across mammals (Cardillo et al., 2005;Ripple et al., 2019).
We also found that the probability of site use for the community of wild species increased with increasing distance from permanent sources of human disturbance, mainly herder houses/settlements, which adds evidence to the general influence of livestock herding on cooccurrence patterns. Notably, such effect was higher than those of other environmental covariates we considered (elevation and slope). On the one hand, this may be the result of contrasting environmental effects at species level compensating each other in the mean effect values, or as species effects are themselves generally weak it may also reflect the relatively uniform habitat that characterises the target landscapes. On the other hand, inter-specific interactions and anthropogenic disturbance might be more important drivers of species distributions than abiotic variables.
While, as predictable, the correlation in occupancy between large-and small-sized livestock was the highest estimated (c = 0.52), it also highlighted differences in site use between the two size classes that support our choice of considering them separately, matching the documented differences in herding practices in the same region, and hence potential impacts on wildlife (Augugliaro et al., 2020). Notably, the snow leopard showed strong evidence for negative correlation in occupancy with both livestock categories, with the stronger being with small-sized livestock. At the same time, temporal overlap in diel activity patterns was higher with large-sized livestock, consistent with them being generally free ranging, while sheep and goats are held in corrals at night (Augugliaro et al., 2020). The negative correlation in occupancy between Siberian ibex and large-sized livestock (c = -0.21) may indicate a detrimental effect of pastoralism on ibex as a consequence of direct competition for grazing Ripple et al., 2014).
Additionally, spatial patterns of interaction among large-sized domestic species and this wild ungulate are possibly sharpened by the high degree of overlap in daily activity patterns (Δ=0.85).
Ibex correlation with small-sized livestock was also negative, though statistically weak, while the positive effect of terrain slope on its site use supports the preference of ibex for rugged habitat (Han et al., 2021).
The positive correlation in occupancy between wolf and large-sized livestock (c =0.22) suggests a role of this carnivore in livestock depredations in the region, especially if considering that such spatial association is not compensated by temporal avoidance (Δ=0.69). This is in accordance with an earlier study in the same area that carried out a focal investigation on predator and prey co-occurrence (Salvatori et al., 2021). We acknowledge that a high degree of spatiotemporal overlap does not necessarily indicate prey preference, but it suggests potential for high encounter rates between carnivores and their prey, which is a key component of prey preference (Fortin et al., 2015;Allen et al., 2021). Conversely, no statistically strong spatial correlation in occupancy was detected for the wolf with herds of small-sized livestock. Overall, results are consistent with Augugliaro et al. (2020), who found wolves to hunt preferably large-rather than small-sized livestock based on interviews to herders.
Marmots and ground squirrels showed positive correlations with both livestock categories, consistent with previous research reporting that these species have likely benefited from the recent increase of livestock in Mongolia (Gankhuyag et al., 2021). However, this apparent tolerance of burrowing mammals towards livestock is at odds with the documented severe decline of Siberian marmots (Marmota sibirica) in Mongolia in the last decades, possibly related to poaching for the fur market (Batbold, 2002;Kolesnikov et al., 2009) and thus requires further investigation.
As for predator-prey interactions, our results indicated a positive spatial correlation of the snow leopard with its primary wild prey in the region, the Siberian ibex .
Similarly, snow leopard site use was positively, although weakly, correlated with that of marmots, which also represent a significant part of its diet (Hacker et al., 2021;Lukarevskiy et al., 2019).
Overall, our results are of relevance for this IUCN-Vulnerable top predator for at least three reasons: (1) they provide an indication of spatial displacement of snow leopards by livestock, especially small-sized, which is concerning given the growth of the cashmere market (Berger et al., 2013;e.g. Tumursukh et al., 2016). (2) Siberian ibex appears as a key driver of occurrence for this carnivore, while there was no strong evidence of an effect of any abiotic covariates on its site use.
(3) Snow leopards in western Mongolia may prey primarily on wild prey and kill livestock opportunistically, even where the abundance of domestic animals is at least one order of magnitude higher (Augugliaro et al., 2020;Sharma et al., 2015;Johansson et al., 2015;Hacker et al., 2021).
These dynamics markedly differ from what was found in other parts of this predator range (Bagchi & Mishra, 2005;Aryal et al., 2014), supporting the notion that the magnitude of human-snow leopard conflicts is highly context-specific (Snow Leopard Network, 2014). No association emerged instead between wolf and ibex, possibly matching the typical preference of these cursorial predators for structurally less rugged habitat than the one selected by ibex (Suryawanshi et al., 2013). The strong evidence of a positive correlation in occupancy between foxes and marmots (c =0.27) could mirror both direct predation by the fox on marmots and the use of marmots' galleries as refuges/resting sites by foxes (Murdoch et al., 2016). Similarly, Pallas's cat displayed a remarkably high spatial association with marmots (c =0.43), whose cavities the cats depend on as denning and resting sites (Ross et al., 2010b), and also a positive association with pikas (c =0.31), considered the small felid's preferred prey (Ross et al., 2010a).
Concerning intra-guild interactions, spatio-temporal overlap between red fox and wolf (c =0.29, Δ=0.63) has been documented in other ecological contexts, where wolves attract foxes through increased opportunities to scavenge (Ferretti et al., 2021), a facilitation also observed in relation to other carnivores such as the Eurasian lynx and the snow leopard (Krofel et al., 2019. As a complementary result supporting co-occurrence, both canids showed similar environmental and anthropogenic covariate effects on site use. On the contrary, no spatial association or avoidance was observed between the two felids, the snow leopard and the Pallas's cat, probably reflecting low dietary niche overlap and infrequent interactions among them. The clear spatio-temporal association between ground squirrels and marmots (c =0.40, Δ=0.91) is not surprising given that they belong to the ecosystem engineering group of colonial-living burrowing animals in grasslands, and they likely present very similar habitat preferences.
As general limitations of our study we acknowledge that: (1) the sampling design was aimed at maximising snow leopard detections, potentially biasing detections of other species and limiting the spatial extent of our inference to the sites sampled; (2) we interpreted occupancy as the probability of site use by the species, given the marked differences among species in home range and movement patterns relative to our sampling design (Neilson et al., 2018), and considering that 'true' occupancy estimates in continuous habitat depend upon the product between the density of the target species and its home range area (Efford & Dawson, 2012); (3) Given the generally high sensitivity of the camera traps we used, the probability of detecting a species was likely largely dictated by its availability for detection (Kéry & Schmidt, 2008), namely by how often a species passed through the camera's detection zone, given that the camera was located within the home range of one or more individuals of that species; (4) the residual correlations in occupancy we estimated may not entirely reflect actual interactions, as they could also be caused by unaccounted differences in habitat use due to missing environmental covariates (Tobler et al., 2019). While limitation 1 is unlikely to affect our inference, given that the sites chosen were used by multiple species of the target community, limitation 2 should be taken into account when interpreting our results. Our sampling design potentially led to overestimate the proportion of area used by species with large home ranges and high mobility (for example snow leopard, wolf and livestock), while underestimating the proportion of area used by species with small home ranges and short distance movements (such as rodents and lagomorphs; Neilson et al., 2018;Efford & Dawson, 2012). This in turn may inflate the positive co-occurrence between species with large home ranges, since they have a higher probability of being detected at the same sampling sites though the home range overlap might, in reality, be limited. Limitation 3 may bear minimal effects on our results, since the duration of our surveys likely allowed to gain a good representation of the asymptotic proportion of area used by each species, i.e., individuals had enough time to cover a consistent proportion of their home range. Regarding limitation 4, we cannot rule out that the inclusion of other abiotic covariates in our model might have changed the output for the biotic relationships. However, the target ecosystem is characterised by low vegetation diversity and most environmental variability derives from topography, which we accounted for by including elevation and terrain slope in our model.
In conclusion, systematic camera trapping proved a valuable tool to study the inter-specific relationships within a community of wild and domestic mammals inhabiting remote and hardly accessible mountain landscapes. In line with our predictions, livestock encroaching into protected areas influenced the whole community of wild mammals, with their average probability of occurrence increasing away from herders' settlements, and a predominant number of potential spatial interactions involving livestock. The sign of the co-occurrence of livestock and wildlife was species-dependent, with larger-sized species that generally showed a tendency to avoid livestock (with the remarkable exception of the wolf that often prey on it), while rodents and lagomorphs tended to positively co-occur with it. As predicted, predators and their focal prey generally had positive correlation in occupancy, while, contrary to our expectations, a number of intra-guild species pairs also showed positive co-occurrence, with no evidence of spatial or temporal niche partitioning, and further, focal investigations are needed to clarify if this pattern is caused by the use of similar resources by species of the same guild. Though pastoralism has co-existed with wildlife for millennia in central Asian grasslands, our findings suggest that the recent dramatic increase in livestock numbers calls for new policies and strategies to decrease the pressure of livestock husbandry on wildlife communities, with special attention on large-sized wild mammals, namely snow leopard and Siberian ibex, which seem to be potentially more sensitive than smaller mammals to the pervasive presence of domestic animals.