Introduction

The conflict derived from livestock predation by carnivores has received much attention because human retaliation has become one of the main threats for carnivore conservation around the world (^{Treves and Karanth 2003}; ^{Inskip and Zimmermann 2009}; ^{Treves and Bruskotter 2014}; ^{Khorozyan et al. 2015}; ^{van Eeden et al. 2017}). Several approaches have been proposed, and various recommendations discussed to reduce the impacts of carnivores (i. e., ^{Graham et al. 2005}; Treves et al. 2006; ^{Kissling et al. 2009}; ^{Pacheco et al. 2008}; ^{Packer et al. 2009}; ^{Hoogesteijn and Hoogesteijn 2011}; ^{Zarco-González et al. 2012}; ^{Miller 2015}). The main global lesson appears to be that, an in-depth assessment of the particularities of the conflict at the local level may be the most promising tool as the basis for an adequate mitigation strategy of any wildlife-human conflict (^{Zapata et al. 2012}; ^{Redpath et al. 2013}).

Here, we report on an analysis of the socio-ecological variables associated to predation by pumas (Puma concolor), and culpeo foxes (Lycalopex culpaeus) on livestock within Sajama National Park (SNP), a protected area located in the Altiplano of Bolivia. The first assessment of the conflict between carnivores and livestock at SNP carried out in 1998 estimated an annual loss of 79 llamas (Lama glama) and 117 alpacas (Lama (Vicugna) pacos; ^{Ribera-Arismendi 1999}). Another research in the same area estimated that the number of llamas and alpacas, as a fraction of the total number of large animals required to maintain a puma population of ~10 individuals at SNP, coincided very closely to the numbers reported by local ranchers in terms of kills by pumas (^{Pacheco et al. 2004}; ^{Gallardo et al. 2009}). SNP is the only protected area in the center of the Bolivian Altiplano that harbors a population of pumas, which is probably connected only to a population in the Chilean Lauca National Park, and a few surrounding areas in Bolivia with low human population density. The main areas where predation occurs at SNP have been mapped along with local herders (^{Ribera-Arismendi 1999}), and apparently, predation was not the major cause of livestock mortality there (^{Zacari and Pacheco 2005}). Therefore, a research assessment of the environmental and livestock management factors associated with the risk of predation is imperative to develop science-based management strategies that balance carnivore conservation with conflict mitigation (^{Zacari and Pacheco 2005}; ^{Weber and Rabinowitz 1996}; ^{Kolowski and Holekamp 2006}; ^{Baker et al. 2008}; ^{Castaño-Uribe et al. 2016}).

A number of socio-ecological variables affect the probability of a carnivore attacking livestock. For example, predation has been shown to increase with livestock density, reductions of native prey, previous attacks in the same area, and distance from human settlements, whereas it tends to decrease closer to roads, and other areas with higher human activity (^{Newmark et al. 1994}; ^{Linnell et al. 2001}; ^{Conover 2002}; ^{Mazzolli et al. 2002}; ^{Novaro et al. 2004}; ^{Woodrofe et al. 2005}; ^{Holmern et al. 2007}; ^{Karlsson and Johansson 2010}; ^{Khorozyan et al. 2015}; ^{Ohrens et al. 2015}). Thus, habitat characteristics are important predictors for predation risk, which in turn depends on the specific predator (^{Stahl et al. 2001}; ^{Soto-Shoender and Giuliano 2011}; ^{Abade et al. 2014}; ^{Zanin et al. 2015}). In particular, for culpeo foxes, livestock predation is significantly affected by both, livestock and native prey densities (^{Pia et al. 2003}). In contrast, for livestock predation by pumas the type of habitat, prey densities (both, native and livestock), and livestock management practices are important factors (^{Polisar et al. 2003}).

Most human communities across the Bolivian Altiplano, including the SNP, have commonly managed the livestock predation conflict by simply killing carnivores. As a consequence, several populations have been eradicated, or severely depleted, including carrion eaters such as the Andean condor (Vultur gryphus), with marked ecological consequences at the community and ecosystem levels, as reported by similar functional losses in other study sites (^{Mech 1998}; ^{Johnson et al. 2001}; ^{Blejwas et al. 2002}; ^{Treves et al. 2002}; ^{Athreya et al. 2004}). Although ^{Gallardo et al. (2009}) estimated that puma population at SNP was around 10 individuals (~1 puma/100 km²), retaliation hunting due to livestock depredation was common in the area (GG, JEJ, pers. obs.). There is no local abundance estimate for culpeo foxes, but we infer that their numbers may be lower than in other regions of their distribution, given the high probability of retaliation hunting by locals and the large home range size (35.9 km²) of a radio-tracked female for ~1 year at SNP (^{Olarte et al. 2009}). This was the second largest estimate for this species, only below the 896 km² home range reported by ^{Jiménez et al. (2001}) and much larger than the 13.9 km² maximum home range size reported for females in Fray Jorge Reserve (^{Salvatori et al. 1999}), both in northern Chile. Therefore, SNP is under a scenario of intense conflict between carnivores and livestock, which is worrisome giving that this protected area harbors protected populations of two of the largest carnivores occurring in the highlands of South America. Consequently, a much-needed contribution to managing the problem entails an analysis of environmental factors that may influence livestock attack. Specifically, in this study we assessed the relative contribution of socio-ecological variables influencing the probability of an attack at the ranch level within SNP. We predicted that the probability of a carnivore attack on a ranch would increase with a) higher livestock densities, b) lower abundance of wild prey, c) good habitat characteristics for carnivores, such as dense vegetation cover, and a rocky structure, and d) low vigilance by ranchers and their watchdogs. Likewise, we predicted that fewer attacks would be associated with e) smaller ranch areas and f) livestock foraging in areas closer to the nearest human settlement. To examine the relative importance of livestock deaths caused directly by carnivores, we also considered new information on other non-predation related factors affecting livestock mortality at the SNP.

Materials and methods

Study area. The Sajama National Park is located on the western portion of the Bolivian Altiplano (-17° 55´ to -18° 15´ S, -68° 41’ to -69° 10’ W), encompassing 100,230 ha at elevations ranging from 4,200 to 6,540 masl. It receives about 400 mm of annual rain, mainly between November and April. Mean annual temperature is 3.4 °C, with a maximum mean of 17.2 °C and a minimum mean of -12.9 °C, during the wet and dry season, respectively. The region is part of the Southern- and the Desert Puna ecoregions of Bolivia (^{Ibisch et al. 2003}), where three vegetation types are distinguished: 1) Keñual-Roquedal, a rocky landscape with sparse forest cover of Polylepis tarapacana trees (locally known as keñua), which grows up to 5,200 m in elevation. Rocks are commonly covered by plants such as Azorella compacta, Parastrephia quadrangularis and Poa asperiflora. This type of habitat is found at higher elevations, where pumas and vicuñas (Vicugna vicugna), the main prey for pumas in terms of biomass (^{Pacheco et al. 2004}), commonly find shelter; 2) Pajonal-Tholar, a steppe type of habitat with softer soil and dominated by Festuca orthophylla, Stipa ichu, and Calamagrostis spp., interspersed with shrubs such as Parastrephia lepidophylla and Baccharis santelicis (locally known as tholas); and 3) Bofedal, a lowland wetland ecosystem dominated by an hydrophilous cushion-type vegetation, dominated by Distichia muscoides and Deyeuxia curvula (^{Ribera-Arismendi 1999}).

Human population at SNP is distributed in five local communities: Caripe, Lagunas, Manasaya, Papelpampa, and Sajama (Figure 1), with a total population of < 800 people estimated a few years before this study (^{Espinoza-Terán 2001}). The main economic activity at SNP is livestock husbandry, mainly of camelids, which are raised in an extensive fashion, encompassing mostly Bofedales, and Pajonal-Tholar, as it is usually the case in the Bolivian Altiplano (^{Cardozo 1985}). At the time of this study, the total livestock at SNP was composed of ~24,000 llamas, ~21,000 alpacas, and ~4,000 sheep distributed in 103 ranches and well over the carrying capacity of the environment (^{Espinoza-Terán 2001}; local ranchers, pers. comm. 2006). As wild camelids, ~3.810 vicuñas roamed free at SNP (^{Asociación Regional de Comunidades Manejadoras de la Vicuña del Parque Nacional Sajama 2006}; ^{Gallardo et al. 2010}).

Figure 1 Study area and spatial distribution of the number of a) puma and b) culpeo fox attacks to livestock on 33 ranches in Sajama National Park, Bolivia, during 2006. 

Three of the four camelid species occurring in Bolivia (^{Anderson 1997}) inhabit SNP. According to local people, the guanaco (L. guanicoe) was eliminated a few decades ago due to hunting pressures. Llamas and alpacas are reared separately, based on species and sex. Females are usually maintained near rural houses at night (but not always in enclosures), and taken to the best pastures (bofedales) during the day, usually below 4,400 masl; although they are sometimes left alone for up to a week. Sheep are usually kept along with female camelids but are confined close to houses at night. Male llamas are maintained away from females, usually in the mountain slopes closer to Keñual-Roquedal (> 4,400 masl). During mating, flocks are joined, and herders closely follow males as they tend to be aggressive toward females and may harm them during courtship. Ranchers use a seasonal movement system for livestock management. During the warm and rainy season, animals are moved to higher areas, and during the cold and dry season, they are taken back to lower elevations (~4,200 masl) to Bofedales. The vicuña, a wild camelid, is the main large prey of pumas at SNP (^{Pacheco et al. 2004}). Vicuñas have been subjected to a live-shearing harvest program for about 18 years.

Livestock husbandry at Sajama National Park. We followed all husbandry activities from January through December 2006 within 33 randomly selected ranches (we did not distinguish between ranches with or without previous carnivore attacks). At each ranch, we interviewed one adult person (either the owner, or the responsible for livestock operations) to obtain the following husbandry information: 1) Seasonal (dry vs wet) spatial distribution of llamas, alpacas, and sheep within the ranch and according to the three vegetation types. 2) The number of people simultaneously looking out for livestock per day. 3) Effort allocated to livestock vigilance as hours/person/week/area, assessed seasonally. 4) Number of shepherd dogs per ranch. 5) Number of domestic animals killed by culpeo foxes, pumas, and other causes (poor-nutrition, diseases, and accidents) during 2006. This information was validated by field confirmation of dead animals.

Prey availability. We estimated abundance of livestock (llamas, alpacas and sheep), vicuñas, Mountain vizcachas (Lagidium viscacia), and small mammals. We standardized data as prey biomass per ranch (kg/area) and calculated prey availability for each ranch as follows: a) Livestock abundance and age structure based on rancher’s knowledge. b) Abundance and biomass for vicuñas according to the SNP annual census carried out in August 2006 (^{Asociación Regional de Comunidades Manejadoras de la Vicuña del Parque Nacional Sajama 2006}). Census procedures include mapping each vicuña group by experienced observers and counting the number of individuals. These locations were digitalized using a GIS to estimate an approximate number of vicuñas at each ranch. c) We estimated the biomass of Mountain vizcacha by sampling six Mountain vizcacha colonies of various sizes, in which we collected fresh feces along a 1-km long and 1-m wide transect per colony. We used the abundance of feces as an index of vizcacha relative abundance, as have been found elsewhere for Mountain vizcachas (^{Walker et al. 2000}) and other Chinchillidae (J. E. Jiménez, unpubl. data). The biomass of Mountain vizcacha was based on abundance indices (fresh fecal pellets/area) obtained in the field and vizcacha density estimates obtained by the method of direct counts, as calibrated for the SNP by ^{Gallardo et al. (2004}). We projected these estimates to the area covered by the Keñual-Roquedal habitat at each ranch, which is the habitat type almost exclusively used by vizcachas at SNP. d) We estimated the density of small mammals (all species combined) at each of the three habitat types by sampling a grid of 10 x 10 Sherman traps, separated 10 m from each other and operated during five nights per habitat, only during the dry season. We baited the traps with a mixture of oats, peanuts, vanilla, and margarine. Traps were checked every morning between 7:00 and 10:30 hrs. Captured animals (only rodents) were marked with a numbered tag for subsequent identification. Abundance was estimated via capture-recapture methods (^{White et al. 1982}).

Estimating probability and frequency of attacks. Several visits to ranches were carried out to confirm all predation events during 2006. For each predation event we recorded: date, prey species, age class of killed individual (young = individuals < 1 year and adults = individuals > 1 year), predator species (according to characteristics described in ^{Shaw et al. 1987}, the experiences of the ranchers, and of our own), exact location, type of habitat where the attack occurred, estimated value of the killed animal, and the owners name. We referred to an attack as a kill by a predator.

Data analyses and modeling. To calculate the biomass of the different animals considered in the study, we used a body mass estimate for each species to extrapolate it at the ranch level, considering the species-specific (except for rodents that were treated as a group) estimated abundances (see above). Thus, for livestock, we used the average body mass estimates provided by local ranchers: 60 kg for adult llamas, 23 kg for young llamas, 45 kg for adult alpacas, 19 kg for young alpacas, and 25 kg for all sheep. Note that these values for llamas and alpacas were very conservative, as these were below the lowest weight for the size ranges estimated for the country: 65-110 kg for llamas, and 50-70 kg for alpacas (^{Cardozo 1985}). The low weights provided by the ranchers were likely the result of the oversized stock in the area. Following the rancher’s opinion on the economic value of killed animals, we used an estimate of $ 13.0 USD for a sheep, $ 48.5 ± 16.9 USD for an adult camelid, and $ 8.6 ± 11.1 USD for a young camelid. These economic values, however, were not used in the modeling analyses. For vicuñas, we used an average body mass of 40 kg for an adult, and 10 kg for a young one (according to unpublished data from ^{Asociación Regional de Comunidades Manejadoras de la Vicuña del Parque Nacional Sajama 2006}). We used an average weight of 1.5 kg for Mountain vizcachas (^{Yensen and Tarifa 1993}). For small mammals, according to our average field estimates, we used a density of 7.3 ind/ha, with an average 30.8 g/ind to estimate biomass for Keñual-Roquedal, and 15.5 ind/ha, with an average of 22.8 g/ind for the biomass estimate in Pajonal-Tholar. We had zero captures in the Bofedal, so this habitat type was excluded from further analysis.

We extrapolated the above estimates to the area covered by each habitat type using GIS tools. We obtained GPS locations (with the owner’s supervision) and plotted each ranch on a vegetation map of SNP (^{Resnikowski 1997}) using ArcView 3.2 (Environmental Systems Research Institute, Redlands, CA, U.S.A.). We then calculated the area for each ranch, and the area covered by each vegetation type within that ranch. We also estimated the distance from the mid-point of a ranch to both, the closest human settlement, and the main small dirt road.

We used the data estimated above to model independently the probability and frequency of an attack on livestock by puma and culpeo fox at the ranch level. We used general linear models (GLMs) and a multimodel inference approach (^{Burnham and Anderson 2002}, 2004) for modeling the likelihood of an attack as explained by the estimated variables, using each ranch as a replicate. Before any general linear modeling was run, however, we conducted pairwise Spearman correlations among all predictor variables to avoid collinearity in the models. All significant correlations with an r _s > 0.9 were considered highly correlated and we eliminated the least relevant variable. In all GLM models, we finally included eight predictor variables: a) ranch area, b) area of the Keñual-Roquedal habitat, c) adult livestock biomass, d) wildlife biomass (vicuña + Mountain vizcacha), e) small mammal biomass, f) distance to the closest human settlement, g) level of husbandry during the dry, and h) during the wet season. Separate estimates of vicuña and Mountain vizcacha biomass were highly correlated with wildlife biomass, so they were excluded from the GLM models. For both puma and culpeo fox we modeled two sets of GLMs based on the type of response variable considered. The first set estimated the probability of an attack using a GLM with binomial distribution and a Logit link function (logistic regression sensu^{Quinn and Keough 2002}). The response variable considered whether a ranch was either attacked (1) or not (0). The second set of models estimated the number of attacks using a GLM with Poisson distribution and Log link function (Poisson regression sensu^{Quinn and Keough 2002}). Here, the response variable was the number of attacks at the ranch level.

Given the large number of models needed to be compared, model selection was based on an autonomous selection process within a multimodel inference approach using the R package glmulti. This approach ranks models from best to worse based on the Akaike Information Criterion corrected for small sample sizes (AICc, ^{Burnham and Anderson 2002}, 200). When more than one candidate model was selected (e. g., ∆AIC ≤ 2) we used model averaging to estimate parameter values proportional to their Akaike weights (w_i). We also calculated a pseudo-R² for the candidate model with the highest number of variables to estimate the proportion of the variance explained. Finally, we compared the livestock husbandry management between seasons with a Wilcoxon test. All statistical analyses were run using SPSS 14.0 (^{SPSS 2005}) and R 3.6.3 (^{R Development Core Team 2020}). Descriptive results are presented as means and one standard deviation (mean ±1 SD).

Results

Characteristics of livestock husbandry. The total number of domestic animals under study was 11,490 for the 33 ranches, including a 50.2 % of llamas, 41.9 % of alpacas, and 7.9 % of sheep, representing ~24 % of the livestock and 34 % of the ranches present at SNP. Mean number of individuals per ranch was 348 ± 233, with an average density of 1.1 ± 1.5 ind/ha (range 104 to 1020 ind/ranch). Mean ranch area was 684.4 ± 744.1 ha. Most of the ranches (~70 %) harbored between 104 to 400 animals, while only 6.1 % of the ranches had more than 800 animals.

At SNP only one family member cares for grazing animals, frequently it was a child. During the dry season, people allocated a smaller amount of time to female livestock vigilance (25.3 ± 12.5 h/week) than during the rainy season (46.3 ± 21.2 h/week, W = 785.1, P < 0.001). Only 21 % of the ranchers had male livestock and this was visited once a week or every other week during the dry season and usually every other day during the rainy season. Ranchers argued that raising male livestock is usually avoided because of recurrent attacks by pumas and because they prefer to allocate more time to other economic activities, such as trading, public transportation, and tourism.

Only 22 (67 %) of the ranches raised dogs, totaling 32 dogs of which only one was a shepherd dog. All ranches had enclosures built near human settlements with a variety of materials including rocks, mud, plants, and Polylepis branches. Enclosures were of different sizes, but 76 % were around 8 m in diameter, and 1.6 m in high (none was higher than 2 m). Most were used mainly for vaccination campaigns during the dry season. Enclosures are not used during the rainy season, because they remain wet and increase the transmission of parasites, and diseases. Sanitary care at SNP was deficient. Only 66 % of the ranchers carried out sanitary campaigns, but those were only partial, as vaccinations were applied only to mothers and youngsters, mostly using inadequate dosages.

Livestock predation by puma and culpeo fox and its correlates. To quantify the number of attacked animals we only considered those that we were able to confirm in situ as killed by a carnivore (84.4 % of all the reported cases). Pumas killed 27 livestock heads, mostly adult camelids (74.1 %), particularly llamas (51.8 %). We did not record any puma attacks on sheep. Of the 156 animals killed by foxes, 90.4 % were young camelids of which 51.6 % were alpacas. Only 18 young sheep were killed by culpeo foxes during our year of monitoring. However, foxes also killed six adult (two years old) llamas and nine adult alpacas during the study (Table 1).

Most attacks by puma (81.5 %) occurred at elevations above 4,500 masl, and away (2,974 ± 3,014 masl) from human settlements; 78 % of these occurred in Keñual-Roquedal habitat type. On the contrary, 89 % of culpeo fox attacks occurred within Pajonal-Tholar, below 4,500 masl, and relatively close (98 ± 78 m) to enclosures. Ranch area was positively correlated with the area of Keñual-Roquedal (r _s = 0.450, d. f. = 32, P = 0.009), distance from human settlements (r _s = 0.741, d. f. = 32, P < 0.001), and livestock abundance (r _s = 0.601, d. f. = 32, P < 0.001).

Livestock losses due to predators. We recorded 183 domestic camelids killed by wild predators during the year of the study (1.6 % of the total heads under study); unleashed dogs killed only one young llama. Foxes were responsible of 85.2 % of the attacks; the rest occurred by pumas. Estimated economic damage totaled ~$ 4,215 USD; ~$ 2,898 USD by foxes, ~$ 1,310 USD by pumas, and ~$ 6.7 USD by dogs. These figures represent the values for 2006 not adjusted for inflation. Predators killed an average of 2.3 ± 0.9 % of the ranches’ heads, with a mean of $127.7 ± 113.2 USD of economic damage per ranch. Around 51 % of the ranches lost < $ 100 USD to predators during the year of study.

Fifty eight percent of the ranches reported predation only by foxes, 12 % only by puma, and 18 % by both predators, with 12 % reporting no losses to predators. The highest number of puma attacks (11 or 40.7 %) was recorded in ranches near the mountains to the NW of SNP, which were also the largest ranches; whereas most attacks by culpeo foxes (44.2 %) occurred in the plains to the S of SNP. Of all the attacks, 84 % occurred during the wet season (January through March; Figure 1). It is interesting to note that ranches with more puma attacks suffered less attacks by foxes and vice versa. An exploratory analysis showed, however, a negative yet marginally significant correlation between the number of attacks by pumas and the number of attacks by culpeo foxes at ranches (r _s = -0.33, n = 33, P < 0.06).

Livestock losses due to other causes. We recorded 353 dead animals that did not show signs of predation, almost two-fold of those killed by predators (Table 1). Among the attributed causes of those deaths were diseases, malnutrition, abortions, and broken limbs. Most deaths were of adult animals (201 vs 152 juveniles), which amounted an estimate of $7,545 USD. This shows that ranches under study lost more animals to other causes (10.7 ± 8.8 animals/ranch) than to both predators combined (5.5 ± 4.3; t = -3.25, n = 33, P = 0.002).

Probability of puma attacks. Of the 256 possible variable combinations (models), the autonomous model selection identified three candidate models (∆AICc ≤ 2) that best explained the probability of a puma attack on ranches (Table 2). These models, in decreasing order of importance, based on averaged parameter values and Log Odd ratios, included ranch area, livestock biomass, and husbandry during the dry season (Table 2). These variables explained 44 % of the variability in the attacks by puma. As the ranch area increased, so did the probability of having a puma attack, being higher than 50 % in ranches larger than 800 ha and almost certain in those having > 1,500 ha (Figure 2 a, b). Likewise, as livestock biomass increased on a ranch, so did the probability of a puma attack, being over 50 % in ranches having ≥ 90 kg/ha of livestock (Figure 2 c, d). Although livestock husbandry during the dry season negatively affected the probability of suffering a puma attack, its effect was rather weak. There was only a 20 % probability of having a puma attack during the dry season in ranches with no husbandry (Figure 2 e, f).

Figure 2 Changes in the probability of puma attacks in relation to socio-ecological variables at Sajama National Park, Bolivia, during 2006. Probabilities are based on odds and log odds ratios from logistic regression GLMs in combination with autonomous multimodel selection. a) Variation in the probability of an attack and b) log odds of an attack due to changes in ranch area (ha). c) Variation in the probability of an attack. d) Log odds of an attack due to changes in livestock biomass (kg/ha) across ranches. e) Variation in the probability of an attack. f) Log odds of an attack due to investment in husbandry during the dry season (h/person/week/area). Grey areas represent the 95 % CI from the model for each variable; dashed lines represent the 50/50 % chance of being attacked; dots are values at the ranch level and thin lines along the top and bottom of each graph show the distribution of values for the predictor variables. 

Frequency of puma attacks. Six candidate models (∆AICc ≤ 2) best explained the number of puma attacks on ranches (Table 2). These models included the density of small mammals, extent of the Keñual-Roquedal habitat, ranch area, and distance to roads as predictors. Although, these four variables explained 64 % of the variability in the frequency of puma attacks, ranch area was the most important variable (95 % CI did not include zero; Table 2). Ranch area was included in all six models. As ranch area increased, the number of puma attacks increased slightly, at a rate of one attack per every 1,000 ha of ranch increase (Figure 3 a,b). In general, the models indicated that once a puma attacked a ranch, the frequency of attacks depended very little on other socio-ecological variables.

Figure 3 Changes in the frequency of puma attacks in relation to socio-ecological variables at Sajama National Park, Bolivia, during 2006. The number and log number of attacks in the model came from Poisson regression GLMs in combination with autonomous multimodel selection. a) Variation in the number of attacks and b) Log of number of attacks (linearized relation) with changes in ranch area (ha). Shaded areas represent the 95 % CI from the model for each variable; dashed lines represent the mean number of attacks and dots are values at the ranch level. 

Probability of culpeo fox attacks. Ten models (∆AICc ≤ 2) best explained the probability of a fox attack on ranches (Table 3). These models included husbandry, both during the wet and during the dry season, livestock biomass, small mammal biomass, vicuña biomass, and the area of Bofedal habitat. Model averaging indicated that husbandry during the wet season was the most important variable, followed by husbandry during the dry season, and the area of Bofedal habitat. The other three variables had odd ratios less than one, indicating that the probability of an attack due to these variables was lower than the probability of not having an attack on a ranch (Table 3). Altogether, the six variables explained only 15 % of the variation. Usually, as husbandry increased in both seasons, the probability of an attack by foxes decreased. Attack probabilities lower than 50 % were reached when people allocated at least from 50 to 60 h/week of vigilance (Figure 4 a-d). Likewise, as Bofedal area increased on a ranch the probability of a fox attack decreased, although the probabilities remained high (> 50 %) regardless of the area (Figure 4 e, f).

Frequency of culpeo fox attacks. Two candidate models (∆AICc ≤ 2) that explained the number of fox attacks on ranches were selected (Table 3). These models included in decreasing order of importance, based on averaged parameter values and their 95 % CI, ranch area, vicuña biomass, livestock biomass and husbandry during the wet season. The effect of husbandry during the wet season was negligible (95 % CI of the parameter estimate included zero; Table 3). Although, the four variables explained 61 % of the variation, their slope estimates had very low values (Figure 5 a-f; Table 3). Nevertheless, as ranch area increased, the number of fox attacks decreased slightly (Figure 5 a, b). The same pattern was found for vicuña- (Figure 5 c, d) and for livestock biomass (Figure 5 e, f). As for puma attacks, our models indicated that once foxes attacked a ranch, the frequency of attacks varied very little with other socio-ecological variables.

Figure 4 Changes in the probability of culpeo fox attacks in relation to socio-ecological variables at Sajama National Park, Bolivia, during 2006. Probabilities are based on odds and log odds ratios from logistic regression GLMs in combination with autonomous multimodel selection. a) Variation in the probability of an attack, and b) log odds of an attack due to changes in Bofedal area (ha). c) Variation in the probability of an attack. d) Log odds of an attack due to changes in livestock biomass (kg/ha) across ranches. e) Variation in the probability of an attack. f) Log odds of an attack due to investment in husbandry during the dry season (h/person/week/area). Shaded areas represent the 95 % CI from the model for each variable; dashed lines represent the 50/50 % chance of being attacked; dots are values at the ranch level and thin lines along the top and bottom of each graph show the distribution of values for the predictor variables. 

Discussion

Most of our predictions turned out as expected, at least for pumas. The probability of a puma attack increased with greater livestock biomass per ranch, yet this probability decreased for culpeo fox attacks. The area covered by dense vegetation and rocky structure (Keñual-Roquedal habitat type) also increased the probability of receiving puma attacks. Increased vigilance by ranchers reduced -though very slightly- the probability of an attack by both predator species. Wild prey abundance, however, affected negatively only the probability of an attack by culpeo foxes.

The importance of a management strategy for predation at SNP becomes clear when considering that the local gross family annual income ranged (at the time of the study) between $290 - 1,500 USD (^{FAO 2005}; C. Espinoza, pers. comm.). This means that the average loss due to predators represented 21.8 ± 19.6 % of the annual family income. This is comparable to what happens in a similar environment in Nepal, where snow leopards (Panthera uncia) killed 2.6 % of sheep holdings, which represented ~25 % of the annual local income (^{Oli et al. 1994}). It seems low, however, when compared to what ^{Valderrama-Vásquez et al. (2016}) reported as an average loss of livestock (sheep and goats) to pumas ($ 110.77 USD per unit of their analysis, that included several ranches) in the Colombian Andes. The high magnitude of the problem at SNP has been pointed out by previous studies (^{Ribera-Arismendi 1999}; ^{Pacheco et al. 2004}; ^{Zacari and Pacheco 2005}). In addition, our mortality data closely resembled those obtained several years ago at SNP by ^{Ribera-Arismendi (1999)}, suggesting that losses due to predation have remained stable. This strengthens our scientific-based recommendations to support the implementation of a management strategy to lessen the conflict (^{van Eeden et al. 2017}) at SNP.

The contrasting effects of the attacks by puma and culpeo fox on livestock can be explained by their presumed foraging behaviors. Pumas, for example, predated mainly adult livestock away from human settlements, in areas toward the mountains, where they usually find shelter; whereas culpeo foxes attacked close to human dwellings, targeted mainly juveniles and included sheep. A likely explanation for this difference is that culpeo foxes remain closer to people to avoid encountering pumas, which commonly prey on them (^{Pacheco et al. 2004}). Note that guardian dogs were uncommon at SNP. Radiotracking studies are needed to address this issue.

Culpeo foxes also preyed on sheep (^{Ribera-Arismendi 1999}), which are usually closer to ranch houses, whereas pumas preyed mainly on vicuñas (^{Pacheco et al. 2004}), which are usually more abundant away from humans at SNP (^{Gallardo et al. 2010}). In more general terms, an interesting interaction between culpeo foxes and pumas may develop when both co-occur. At least in the Bolivian Andes, where pumas are present in areas with livestock, a strong conflict develops due to predation, which is further exacerbated with livestock predation by culpeo foxes (^{Pacheco et al. 2010}). The fact that pumas and culpeo foxes at SNP did not usually use the same habitats for preying upon livestock suggests that although these predators share the same landscapes, they tended to partition habitat types. We hypothesize that this is driven by intra-guild predation of pumas on foxes (^{de Oliveira and Pereira 2014}).

Two other findings are notable from our results: i) dogs do not usually attack livestock at SNP. Dogs may be important predators for wildlife elsewhere in South America (^{Schüttler et al. 2018}; ^{Rodriguez et al. 2019}), and they are usually blamed for attacking livestock in other regions of Bolivia (^{Aliaga-Rossel et al. 2012}). Although dogs are well-known allies in reducing carnivore attacks (^{Gonzales et al. 2012}; ^{Novaro et al. 2017}; ^{Spencer et al. 2020}), Sajama ranchers are not prone to raising dogs, because they are supposedly a source of parasites. Therefore, dogs at SNP are rather rare and problems with unleashed dogs are seldom reported. ii) Losses of domestic camelids due to diseases and undernourishment took at least twice as many heads as those preyed by both predators combined. This is in line with ^{Zacari and Pacheco’s (2005}) findings, who, based on clinical examinations and theoretical mortality rates, estimated that diseases would kill two to six times more livestock than predators at SNP, a figure very close to what was estimated in this study. Furthermore, rural, non-technologically oriented rearing operations for domestic camelids (llamas and alpacas) may usually expect between 45 to 55 % fertility (^{Cardozo 1985}; ^{Quispe et al. 2009}). Mortality at weaning in domestic camelids is quite high though (6 to 8 %) and adult annual mortality may reach up to 10 % (^{Cardozo 1985}) and can be as high as 50 % under drought conditions (^{Quispe et al. 2009}). It is notable that in a region nearby SNP (Turco, Bolivia), droughts have been reported by local people as the most important cause of death for domestic camelids, whereas pumas and culpeo foxes were second and third in importance, respectively, well above the importance of diseases (^{Romero-Muñoz et al. 2016}). Under this scenario, a 1.6 % mortality due to predation at SNP does not seem to be high, at least as compared to losses due to pumas in central Argentinean rangelands (medians: 0.1 to 1.0 % of cattle holdings plus 3.3 to 10.4 for sheep; ^{Guerisoli et al. 2017}) or 0.94 % and 0.39 % of the sheep stock by pumas and culpeos, respectively, in Aysen, southern Chile (J. Jiménez et al. unpublished data).

Figure 5 Changes in the frequency of culpeo fox attacks in relation to socio-ecological variables at Sajama National Park, Bolivia, during 2006. The number and log number of attacks in the model came from Poisson regression GLMs in combination with autonomous multimodel selection. a) Variation in the number of attacks and b) log of number of attacks (linearized relation) with changes in ranch area (ha). c) Variation in the number of attacks. d) Log number of attacks due to changes in livestock biomass (kg/ha) across ranches. e) Variation in the number of attacks. f) Log number of attacks due to changes in vicuña biomass (kg/ha) across ranches. Shaded areas represent the 95 % CI from the model for each variable, dashed lines represent the mean number of attacks and dots are values at the ranch level. 

Probability of attacks by predators. Although intense husbandry may reduce the probability of livestock suffering an attack by puma, ranch area and livestock biomass better predicted that likelihood. Large ranches harbored more livestock (number of animals), which means more potential prey, under less human supervision. Ranches > 1,500 ha in SNP have a probability close to one for receiving an attack from pumas. Furthermore, the number of livestock killed by pumas on ranches larger than 1,500 ha rises exponentially. Finally, larger ranches are usually composed of larger areas of Keñual-Roquedal habitat type, which, according to our long-time, but non-quantitative observations, as well as by locals’ knowledge, is puma’s preferred habitat for shelter at SNP. In other regions, predation by pumas has also been reported to be higher when livestock moves closer to forests and water, away from humans (^{Mazolli et al. 2002}). The preference for pumas for rocky and vegetation cover is likely, in addition to the fact that they usually avoid areas inhabited by humans (^{Dellinger et al. 2020}). Similar patterns of a positive correlation between ranch size and/or livestock numbers and predator attacks have been widely reported for a large number of predator species both in America (^{Mech et al. 2000}; ^{Treves et al. 2004}; ^{Bradley and Pletscher 2005}; ^{Romero-Muñoz et al. 2016}; ^{Sarmiento-Giraldo et al. 2016}), as elsewhere (^{Kaartinen et al. 2009}; ^{Hanley et al. 2018}; ^{Mijiddorj et al. 2018}).

Culpeo foxes seem to respond very differently than pumas to the same socio-ecological variables evaluated. The single most important variable that may be controlled through management practices is husbandry, whose rate of increase is related to a decrease in the probability of receiving an attack from culpeo foxes. Husbandry and herding have been repeatedly found to reduce predator attacks in previous studies elsewhere (^{Conover 2002}; ^{Treves et al. 2002}; ^{Treves and Karanth 2003}; ^{Ogada et al. 2003}; ^{Treves et al. 2006}; ^{Mijiddorj et al. 2018}; ^{Adhikari et al. 2020}). Two other variables that were negatively related to the likelihood of fox attacks were livestock biomass and area covered by Bofedales. Although livestock uses preferentially the Bodefales to forage, it should be noted that our capture success for small mammals at the Bofedales was zero, which may partially explain the absence of culpeo foxes roaming around this type of habitat. In accordance, the study by ^{Olarte et al. (2009}) at SNP suggests that culpeo foxes do not select the use of Bofedales. Furthermore, Bofedales have very low vegetation cover rendering this a highly open and exposed habitat, increasing the visibility of a predator such as a culpeo fox to humans and pumas.

In a review of > 120 estimates of livestock loss to predation by carnivores across five continents, ^{Baker et al. (2008}) noted that it usually ranged between < 1 and 5 % of the total livestock holdings. In 6 % of the cases, however, losses by predation reached > 20 % of the stocks. Compared to other case studies, our yearly loss estimate of 1.6 % of the stock to predation was higher than the maximum of 0.59 % reported in France (^{Stahl et al. 2001}) as a yearly loss of sheep due to Eurasian lynx (Lynx lynx) predation, but lower than the 4.6 % combined losses to jaguars (Panthera onca) and pumas on livestock at four sites in Costa Rica (^{Amit and Jacobson 2017}), and even lower than the 10 % herd size losses due to Ethiopian wolves (Canis simensis, ^{Eshete et al. 2018}) or the 23.9 % losses of livestock due to wolves (Canis lupus) and snow leopards (Panthera uncia) in Mongolia (^{Mijiddorj et al. 2018}). A loss of 0.25 to 5.0 % of the livestock was reported in North American and European studies, whereas it tended to be higher (< 1 to ~10 % of the total stocks) in African sites (^{Thirgood et al. 2005}). Under certain circumstances, however, losses can be as high as 78 % as has been the case for goats lost to pumas in southern Brazil, after carnivores were protected (^{Mazzolli et al. 2002}), or the 40 % of sheep lost to culpeo foxes in some ranches in Argentina, at times when native prey was scarce (^{Novaro et al. 2004}).

The absence of a relationship between native prey biomass and probability of an attack from predators found in our study coincides with the general patterns found by ^{Ugarte et al. (2019}). The probability of suffering an attack by pumas and culpeo foxes reached almost 100 % for large ranches, irrespective of the wild prey biomass available. Therefore, it seems reasonable to conclude that the probability of suffering an attack at any single ranch in the SNP is quite high, but the accumulated losses to predation are rather low, as compared to other sites around the world.

Although the numerical losses at SNP may seem low, their effect on the local economy was not. Losses to predators represented between 8 to 44 % of the annual family income at SNP. Five studies from Africa and Asia reviewed by ^{Suryawanshi et al. (2017}) revealed that up to 50 % of the annual income may be lost to predators. Thus, our figures are not the highest. An important fact to bear in mind, however, is that the proportion of a family’s income that is lost to predators, may be important when one tries to implement a mitigation measure to reduce predation (^{Amit and Jacobson 2017}). In fact, our long-term experiences at SNP indicate that local people have very low credibility on outsiders (scientists among them) when discussing mitigation measures for wildlife depredation. This certainly adds another dimension to diminishing the conflict.

The probability of suffering an attack from pumas increased with ranch size, livestock biomass, and closeness to the rocky and well-vegetated habitats embodying shelter for pumas. Husbandry also reduced the probability of a puma attack, but only marginally. Culpeo fox attacks were more likely in ranches with less livestock biomass and a larger area covered by Bofedal. Greater attention on the herds by humans contributed to reducing the probability of a fox attack.

The most important source of mortality for domestic camelids at SNP was not predation by either pumas or foxes, but by diseases and malnutrition. This fact prompts us to pose, as a hypothesis (see ^{Treves and Santiago-Ávila 2020}), that livestock predation by carnivores at SNP may be acting as compensatory mortality, if weak animals are killed more often (i. e., ^{Bradley and Pletscher 2005}). We do not have data to test this hypothesis yet.

Given that the livestock standing biomass at SNP is higher than in most of the Bolivian Altiplano (^{Cardozo 1985}; ^{MAPZA-GTZ 2001}), increasing health care may only make things worse concerning carrying capacity (already surpassed) and, hence, overall health indices for camelids at SNP.

Although we do not have data from other regions of Bolivia to compare with, our findings on the causes of mortality at SNP let us venture to propose a few recommendations:

a) Standing livestock biomass should be reduced, to increase overall health of all domestic camelids. As an indirect consequence, fewer predation events by pumas are expected and degraded vegetation may improve.

b) Once livestock health indices have improved, the critical and feasible measure to reduce predation seems to be the increase in husbandry practices. Although our modeling shows an overall weak effect of husbandry, we must consider, for instance, that male llama groups are left alone on prime puma habitat and are only visited two to four times a month. One possibility is to establish communal herders, which should be compensated relative to the time invested in tending livestock, while providing livestock owners more free time to focus on other economic activities.

c) The use of shepherd dogs, which is highly recommended by several studies elsewhere, including Argentina (^{Gonzales et al. 2012}; ^{Novaro et al. 2017}) should be an additional strategy to be implemented at SNP. Although local people do not like to raise dogs, their reasoning behind this attitude can be overturned if a program includes the vaccination of shepherd dogs. Clearly, an economic cost-benefit analysis of raising and caring for a shepherd dog (i. e., ^{Saitone and Bruno 2020}) should be carried out before implementing such a program at SNP.

d) Other potential measures that can be implemented are the use of conditioned taste aversion and exploring a predator friendly market for camelid meat from SNP (^{Cowan et al. 2000}; ^{Aquino and Falk 2001}; ^{Pacheco et al. 2008}). Implementing the latter would provide better revenues for herders while accepting some losses to predators.

e) Prohibit hunting of pumas as a measure to reduce attacks, given that the puma population at SNP may be too small to withstand any hunting pressure (^{Gallardo et al. 2010}). If hunting for culpeo foxes is to be considered, a population assessment should be done before any actions can be taken, and learning from experience, such as the one with > 30 years of culpeo control in Argentina (^{Funes et al. 2006}) should be a priority.