DOI: 10.5061/dryad.n5tb2rbv6
title: Spatial separation of prey from livestock facilitates coexistence of a specialized large carnivore with human land use

description: There is an increasing emphasis in conservation strategies for large
 carnivores on facilitating their coexistence with humans. Justification
 for coexistence strategies should be based on a quantitative assessment of
 currently remaining large carnivores in human-dominated landscapes. An
 essential part of a carnivore’s coexistence strategy has to rely on its
 prey. In this research, we studied snow leopards Panthera uncia whose
 habitat mainly comprises human-dominated, unprotected areas, to understand
 how a large carnivore and its primary prey, the bharal Pseudois nayaur,
 could coexist with human land use activities in a large proportion of its
 range. Using a combination of livestock census, camera trapping and
 wildlife surveys, across a broad gradient of livestock grazing intensity
 in a 363,000 km2 landscape on the Tibetan Plateau, we found no evidence of
 livestock grazing impacts on snow leopard habitat use, bharal density and
 spatial distribution, even though livestock density was 13 times higher
 than bharal density. Bharal were found to prefer utilizing more rugged
 habitats at higher elevations with lower grass forage conditions, whereas
 livestock dominated in flat valleys at lower elevations with higher
 productivity, especially during the resource-scarce season. These findings
 suggest that the spatial niche separation between bharal and livestock,
 together with snow leopards’ specialized bharal diet, minimized conflicts
 and allowed snow leopards and bharal to coexist in landscapes dominated by
 livestock grazing. In recent years, reduced hunting and nomadic herder’s
 lifestyle changes towards permanent residence may have further reinforced
 this spatial separation. Our results indicated that, for developing
 conservation strategies for large carnivores, the niche of their prey in
 relation to human land-use is a key variable that needs to be evaluated.
description: # Spatial separation of prey from livestock facilitates coexistence of a
 specialized large carnivore with human land use
 [https://doi.org/10.5061/dryad.n5tb2rbv6](https://doi.org/10.5061/dryad.n5tb2rbv6) This dataset contains all data used in the data analysis procedure of this study, including SL\_watershed0.csv to construct the GLMM for modelling snow leopard habitat use, and BSden\_unit0.csv to construct the GLMM for modelling bharal density. The other three files, BS.occupancy.csv, BS.dis.csv, BS.sea.csv, BS.site.csv were used to construct occupancy model for bharal spatial distribution. As authors, we confirm the data didn't contain any location data for any of the species involved, thus won't affect the conservation of snow leopard. ## Description of the data and file structure Abbreviations in SL\_watershed0.csv: ID (valley ID), SL (snow leopard detection number), CD (camera working day), BS (bharal density/km^2), LS (livestock grazing intensity in sheep unit/km^2),EVI (mean grassland EVI), WP (mean winter precipitation),WT (mean winter temperature), DEM (mean elevation), VRM (mean ruggedness), Slope (mean slope), House (house density/km^2), Rock (rock coverage), Site (study site name) For bharal occupancy model: BS.occupancy.csv: Bharal detection status in each 2 by 2 km grid across 8 surveys. BS.site.csv: occupancy covariates. Grid (grid ID), LS (livestock presence), Rock (rock coverage), DEM (elevation), EVI (grassland EVI), Slope, VRM (ruggedness), Site (study site name) BS.sea.csv and BS.dis.csv are two detection covariates used in occupancy model, one indicating the survey seasona and one indicating the distance from the transects to bharal herds. ## Code used for plotting the figures: library(plyr) library(reshape2) library(ggplot2) library(hrbrthemes) #preparing data setwd("D:/phD/DATA\_ANALYSIS/Raw\_data") bs\_ls&lt;-read.csv("BS\_LS\_loc.csv", sep=",") head(bs\_ls) bs\_ls$Species&lt;-as.factor(bs\_ls$Species) bs\_ls$Season&lt;-as.factor(bs\_ls$Season) bs\_ls$Groups&lt;-as.factor(bs\_ls$Groups) attach(bs\_ls) library(car) library(rgl) install.packages("barplot3d") library(barplot3d) scatter3d\(Elevation\~EVI\+Slope\|Species\, data=bs\_ls\, axis\_ticks=T\, surface=F\) melted &lt;- melt(bs\_ls, id.vars=c("Species")) com&lt;-ddply(melted, c("Species", "variable"), summarise, mean = mean(value), sd = sd(value), sem = sd(value)/sqrt(length(value))) com$upper&lt;-com$mean+1.96\*com$sem com$lower&lt;-com$mean-1.96\*com$sem EVI&lt;-subset(com,variable=="EVI") Slope&lt;-subset(com,variable=="Slope") Elevation&lt;-subset(com,variable=="Elevation") \# This part below was just to play around with the data to see how it looks \# scatter plot with confidence intervals p1&lt;-ggplot(EVI, aes(Species, mean))+xlab("Species")+ylab("EVI")+ geom\_pointrange(aes(ymin = lower, ymax = upper))+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p2&lt;-ggplot(Slope, aes(Species, mean))+xlab("Species")+ylab("Slope")+ geom\_pointrange(aes(ymin = lower, ymax = upper)) + theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p3&lt;-ggplot(Elevation, aes(Species, mean))+xlab("Species")+ylab("Elevation")+ geom\_pointrange(aes(ymin = lower, ymax = upper))+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) require(gridExtra) grid.arrange(p1, p2, p3, ncol=3) \# density plot p1 &lt;- ggplot(data=bs\_ls, aes(x=EVI, group=Species, fill=Species)) + geom\_density(adjust=1.5, alpha=.4) + theme\_ipsum(axis\_title\_size = 12, axis\_title\_face = "plain", grid = F, axis\_col = "#cccccc", axis = T, ticks = T) p2 &lt;- ggplot(data=bs\_ls, aes(x=Slope, group=Species, fill=Species)) + geom\_density(adjust=1.5, alpha=.4) + theme\_ipsum(axis\_title\_size = 12, axis\_title\_face = "plain", grid = F, axis\_col = "#cccccc", axis = T, ticks = T) p3 &lt;- ggplot(data=bs\_ls, aes(x=Elevation, group=Species, fill=Species)) + geom\_density(adjust=1.5, alpha=.4) + theme\_ipsum(axis\_title\_size = 12, axis\_title\_face = "plain", grid = F, axis\_col = "#cccccc", axis = T, ticks = T) grid.arrange(p1, p2, p3, nrow=3) \# Script for figure 2 setwd("D:/Manuscripts/Land sharing in mountainous refuge") c&lt;-read.csv("Coefficients\_no GI.csv", sep=",") head(c) c$Variable &lt;- factor(c$Variable, levels = rev(c("Grazing intensity","Bharal density","Winter temperature","Winter precipitation","Elevation","Vision coverage"))) p1&lt;-ggplot(c, aes(Coefficients, Variable, colour=Model, order=ID))+labs(fill="Models")+ scale\_colour\_discrete(limits=c("Snow leopard habitat use", "Bharal density"))+ xlab("Standardized coefficient with 95% CI")+ylab("Covariates")+ geom\_pointrange(aes(xmin = Lower, xmax = Upper),position=position\_dodge(width=0.5))+ geom\_vline(xintercept = 0, linetype = "dashed")+ xlim(-1.2, 2)+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p1 bs&lt;-read.csv("BSoccu.csv", sep=",") head(bs) bs$X &lt;- factor(bs$X, levels = rev(c("p (Survey season)","p (Distance to transect)","psi (Livestock presense)","psi (Elevation)","psi (Slope)","psi (EVI)"))) p2&lt;-ggplot(bs, aes(Coefficient, X, colour=Model)) + scale\_colour\_manual(values=c("goldenrod4", "goldenrod2"))+xlab("Standardized coefficient with 95% CI")+ylab("Covariates")+ geom\_pointrange(aes(xmin = Lower, xmax = Upper),position=position\_dodge(width=0.5))+ geom\_vline(xintercept = 0, linetype = "dashed")+ theme(panel.grid.major = element\_blank(), panel.grid.minor = element\_blank()) p2 require(gridExtra) grid.arrange(p1, p2, nrow=2) #PCA analysis for niche partitioning visualization \# Figure 3b bs\_ls&lt;-read.csv("ALL\_BS\_LS\_PCA2.csv", sep=",") head(bs\_ls) library(ggfortify) df &lt;- bs\_ls[1:8] pca\_res &lt;- prcomp(df, scale. = TRUE) autoplot(pca\_res) p1 &lt;- autoplot(pca\_res, data = bs\_ls, colour = 'Species',shape= 'Season',loadings = TRUE, loadings.colour = 'blue', loadings.label = TRUE, loadings.label.size = 4) + stat\_ellipse(geom="polygon", aes(fill = Species\_Season), alpha = 0.2, show.legend = FALSE, level = 0.95)+ theme\_minimal() + theme(panel.grid = element\_blank(), panel.border = element\_rect(fill= "transparent")) p1 #PCA with FactoMineR library(FactoMineR) res.pca = PCA(bs\_ls, quali.sup=9:11, graph=T) summary(res.pca) plot.PCA(res.pca, label = c("ind.sup","quali","var"), habillage="Species") require(factoextra) require(ggplot2) require(tidyr) require(dplyr) require(MASS) require(reshape2) require(cowplot) \# extract pc scores for first two component and add to dataframe bs\_ls$pc1 &lt;- res.pca$ind$coord[, 1] # indexing the first column bs\_ls$pc2 &lt;- res.pca$ind$coord[, 2] # indexing the second column bs\_ls$pc3 &lt;- res.pca$ind$coord[, 3] # indexing the third column write.csv(bs\_ls[,11:14], "PC\_BS\_LS.csv") \# At last we didn't use the figures produced here p &lt;- ggplot(data = bs\_ls, aes(x = pc1, y = pc2, color = Species, shape = Season)) + geom\_hline(yintercept = 0, lty = 2) + geom\_vline(xintercept = 0, lty = 2) + geom\_point(alpha = 0.8) p &lt;- p + stat\_ellipse(geom="polygon", aes(fill = Species), alpha = 0.2, show.legend = FALSE, level = 0.95) + xlab("PC 1 (27.67%)") + ylab("PC 2 (17.95%)") + theme\_minimal() + theme(panel.grid = element\_blank(), panel.border = element\_rect(fill= "transparent")) p \# Wilcoxon test on PC1 and PC2 wilcox.test(pc1\~Species, data=bs\_ls) wilcox.test(pc2\~Species, data=bs\_ls) #nicheRover for niche partitioning visualization (Figure 3c, d) library(nicheROVER) bsls&lt;-read.csv("PC\_BS\_LS.csv", sep=",") head(bsls) aggregate(bsls[2:3], bsls[1], mean) \# generate parameter draws from the "default" posteriors of each fish nsamples &lt;- 1e3 system.time({ bsls.par &lt;- tapply(1:nrow(bsls), bsls$Species, function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3])) }) \# 2\-d projections of 10 niche regions clrs &lt;- c("black", "red", "blue", "orange") # colors for each species nsamples &lt;- 10 bsls.par &lt;- tapply(1:nrow(bsls), bsls$Species, function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3])) \# format data for plotting function bsls.data &lt;- tapply(1:nrow(bsls), bsls$Species, function(ii) X = bsls[ii,2:3]) p2 &lt;- niche.plot(niche.par = bsls.par, niche.data = bsls.data, pfrac = .05, iso.names = expression(PC1, PC2), col = clrs, xlab = expression("scores")) \# niche overlap plots for 95% niche region sizes nsamples &lt;- 1000 bsls.par &lt;- tapply(1:nrow(bsls), bsls$Species, function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3])) \# Overlap calculation\. use nsamples = nprob = 10000 \(1e4\) for higher accuracy\. \# the variable over\.stat can be supplied directly to the overlap\.plot function over.stat &lt;- overlap(bsls.par, nreps = nsamples, nprob = 1e3, alpha = c(.95, 0.99)) #The mean overlap metrics calculated across iteratations for both niche #region sizes (alpha = .95 and alpha = .99) can be calculated and displayed in an array. over.mean &lt;- apply(over.stat, c(1:2,4), mean)\*100 round(over.mean, 2) over.cred &lt;- apply(over.stat\*100, c(1:2, 4), quantile, prob = c(.025, .975), na.rm = TRUE) round(over.cred[,,,1]) # display alpha = .95 niche region \# Overlap plot\. Before you run this\, make sure that you have chosen your alpha level\. clrs &lt;- c("black", "red", "blue", "orange") # colors for each species over.stat &lt;- overlap(bsls.par, nreps = nsamples, nprob = 1e3, alpha = .95) p3 &lt;- overlap.plot(over.stat, col = clrs, mean.cred.col = "turquoise", equal.axis = TRUE, xlab = "Overlap Probability (%) -- Niche Region Size: 95%") \# posterior distribution of \(mu\, Sigma\) for each species nsamples &lt;- 1000 fish.par &lt;- tapply(1:nrow(fish), fish$species, function(ii) niw.post(nsamples = nsamples, X = fish[ii,2:4])) \# posterior distribution of niche size by species bsls.size &lt;- sapply(bsls.par, function(spec) { apply(spec$Sigma, 3, niche.size, alpha = .95) }) \# point estimate and standard error rbind(est = colMeans(bsls.size), se = apply(bsls.size, 2, sd)) \# boxplots clrs &lt;- c("black", "red", "blue", "orange") # colors for each species p4 &lt;- boxplot(bsls.size, col = clrs, pch = 16, cex = .5, ylab = "Niche Size", xlab = "Species\_Season") require(gridExtra) grid.arrange(p1, p2, p3, nrow=2, ncol=2)