Spatial separation of prey from livestock facilitates coexistence of a specialized large carnivore with human land use
Data files
Sep 25, 2023 version files 63.71 KB
-
BS.dis.csv
14.05 KB
-
BS.occupancy.csv
6.47 KB
-
BS.sea.csv
6.48 KB
-
BS.site.csv
14.06 KB
-
BSden_unit0.csv
6.47 KB
-
README.md
10.26 KB
-
SL_watershed0.csv
5.91 KB
Abstract
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.
README: Spatial separation of prey from livestock facilitates coexistence of a specialized large carnivore with human land use
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<-read.csv("BS_LS_loc.csv", sep=",")
head(bs_ls)
bs_ls$Species<-as.factor(bs_ls$Species)
bs_ls$Season<-as.factor(bs_ls$Season)
bs_ls$Groups<-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 <- melt(bs_ls, id.vars=c("Species"))
com<-ddply(melted, c("Species", "variable"), summarise,
mean = mean(value), sd = sd(value),
sem = sd(value)/sqrt(length(value)))
com$upper<-com$mean+1.96*com$sem
com$lower<-com$mean-1.96*com$sem
EVI<-subset(com,variable=="EVI")
Slope<-subset(com,variable=="Slope")
Elevation<-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<-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<-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<-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 <- 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 <- 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 <- 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<-read.csv("Coefficients_no GI.csv", sep=",")
head(c)
c$Variable <- factor(c$Variable, levels = rev(c("Grazing intensity","Bharal density","Winter temperature","Winter precipitation","Elevation","Vision coverage")))
p1<-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<-read.csv("BSoccu.csv", sep=",")
head(bs)
bs$X <- factor(bs$X, levels = rev(c("p (Survey season)","p (Distance to transect)","psi (Livestock presense)","psi (Elevation)","psi (Slope)","psi (EVI)")))
p2<-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<-read.csv("ALL_BS_LS_PCA2.csv", sep=",")
head(bs_ls)
library(ggfortify)
df <- bs_ls[1:8]
pca_res <- prcomp(df, scale. = TRUE)
autoplot(pca_res)
p1 <- 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 <- res.pca$ind$coord[, 1] # indexing the first column
bs_ls$pc2 <- res.pca$ind$coord[, 2] # indexing the second column
bs_ls$pc3 <- 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 <- 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 <- 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<-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 <- 1e3
system.time({
bsls.par <- 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 <- c("black", "red", "blue", "orange") # colors for each species
nsamples <- 10
bsls.par <- tapply(1:nrow(bsls), bsls$Species,
function(ii) niw.post(nsamples = nsamples, X = bsls[ii,2:3]))
# format data for plotting function
bsls.data <- tapply(1:nrow(bsls), bsls$Species, function(ii) X = bsls[ii,2:3])
p2 <- 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 <- 1000
bsls.par <- 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 <- 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 <- apply(over.stat, c(1:2,4), mean)*100
round(over.mean, 2)
over.cred <- 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 <- c("black", "red", "blue", "orange") # colors for each species
over.stat <- overlap(bsls.par, nreps = nsamples, nprob = 1e3, alpha = .95)
p3 <- 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 <- 1000
fish.par <- 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 <- 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 <- c("black", "red", "blue", "orange") # colors for each species
p4 <- 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)