Methods
2.1 Determining pika occurrence - Appropriate microhabitats
were visually inspected for signs of pika activity (fecal pellets,
haypiles, and burrows). A 10m*10m quadrat was laid centered around each
identified site. In each quadrat, variables about rocks, vegetation, and
topographic features were documented (Table S1). In addition to each
occurrence plot, a random paired plot was laid 50m away (compass
direction determined by a uniformly generated random number bound
between 1 and 360). For social pikas, colony boundaries were identified
using outermost burrows, and paired plots were laid outside the colony.
2.2 Extraction of high-resolution topographic information from
remotely sensed data - For mapping altitude, a Digital Elevation Model
(30m) for the tiles of our interest was downloaded from the Bhuvan,
Indian Geo Platform of ISRO (https://bhuvan.nrsc.gov.in/home/index.php).
The raster was further processed in ArcMap 10.5 (ESRI 2022) to give
slope and aspect, which represent relevant ecological attributes of a
landscape that animals are known to use to select sites for residence.
NDVI values (30m) for each plot location were obtained for the months of
July-August (Summer) and November-December (Winter) dating back nine
years (2013 to 2021), from Landsat 8 imagery (USGS) with the least cloud
cover, using Google Earth Engine
(Gorelick et al. 2017).
These were further processed to provide additional variables concerning
NDVI (Table S1).
2.3 Sampling sites - Based on logistics and pilot surveys,
sampling locations were chosen across elevational distribution ranges of
the three pika species (Figure 2, Figure 3). Since we aimed to
understand ecological factors that shape pika presence at a broad scale,
plots were laid in different locations for each species (Figure 2)
(South-East Ladakh - 88 plots; North-West Ladakh - 66). Sampling was
done in zones of sympatry and supplemented with sampling in zones of
allopatry to capture the niche breadth of all species. The regions
sampled included South-East Ladakh (part of the Qinghai-Tibetan plateau;
4000kmsq) and North-West Ladakh (Ladakh and Zanskar ranges; 8000kmsq).
2.4 Generalised Linear Models (GLMs) - To identify ecological
factors that support the presence of different species of pikas at a
broader scale, quadrats sampled across species and sites (n=148) were
pooled and scored for a binary presence/absence variable for each
species (for instance presence and absence plots of ON andOL were treated as absence points for OM ). Such an
approach is helpful as it is informative regarding where a species is
likely or unlikely to be present. We modeled the presence of pikas at
sites using logistic regressions with plot variables (Table S1) as
predictors using the ‘glm’ function in the R package ‘arm’
(R core team 2021, Gelman
and Su 2021).
Landscape features and vegetation structure were noticeably different
across sites (Figure 9), with huge-sized boulder rocks at low elevations
in North-West Ladakh (Thanh
et al. 2010), including sea buckthorn along rivers and in South-East
Ladakh, the extensive cover of Caragana , cushion plants, and
small-sized scree rocks at high elevations. We, therefore, modeled the
presence of these species at three different spatial extents: a)
South-East and North-West Ladakh (controlling for differences in
landscape and vegetation), b) Ladakh (across the entire landscape).
Since only OM and ON are found in North-West Ladakh, only the presence
of these species was modeled at this spatial scale. At the South-East
Ladakh scale, the occurrence of all three species was modeled as they
co-occur in this landscape.
Due to problems associated with sample size, all variables additively
could not be subject to a stepwise/backward selection algorithm. To
resolve this problem, we use a modified version of the purposeful
selection framework
(Bursac et al. 2008,
Zhang 2016). The workflow used is depicted below:
- Predictor variables that affected presence were identified using
univariate logistic regressions and plots (one predictor and one
binary response variable). All variables that explained the occurrence
of a species in these models (p<0.05) were additively used
to model occurrence (Model 1: Statistical model). In addition, all
variables with a p-value cut-off of 0.25 and other variables that were
biologically relevant (based on field observations) were additively
used to model the occurrence (Model 2: Purposeful model) following
Zhang (2016). Finally, a
third model was constructed with predictors (p-value < 0.25
in uni-variate regressions) binned into different ecological classes
(Table S2) to depict the different landscape-level features (Model 3:
sub-model) that could influence pika presence.
- A stepwise regression (StepAIC) was run on all three models to pick
models that best explain the occurrence (based on the three species of
interest). When AIC values varied marginally across sub-models in
ecological classes, all predictors across models were considered for
further analysis.
- For the sub-model alone, variables from all the best sub-models across
different ecological classes were additively used to model presence
using a stepwise regression for the second time. While doing so, care
was taken to add/replace biologically relevant variables (based on
field experience and literature) even if they were dropped during the
sub-modeling process.
- Best models across modeling approaches were compared and analyzed.
Finally, biologically relevant interactions (based on field
observations) were evaluated based on AIC scores. Goodness of fit
tests (Hosmer and Lemeshow 1989) were performed to assess model fits.
2.5 Comparison of niche overlaps and niche width between species
at broad and regional scales - A dataset was constructed by filtering
out all absence points and including only variables that were both
biologically (based on natural history observations) and statistically
important (based on univariate regressions for each species) for the
presence of each species (OM - variables related to rocks and geography;
ON - variables related to vegetation and shrubs in particular; OL -
variables related to vegetation, elevation) (variables highlighted in
blue on Table S2). This dataset was subjected to a PCA analysis to
reduce dimensionality and to visualize the variation along eigenvectors.
The PCA revealed the separation of broad niches along different axes of
rocks, vegetation, elevation, and NDVI. Variables related to vegetation
(perc_vegetation, height_tallgrass_reeds, height_shrub, average NDVI
in summer) and rocks (perc_rocks) loaded onto PC1, while PC2 and PC3
were driven byloaded on by only one variable, elevation and variance in
winter NDVI respectively (cos2 >0.5).
Since three PC axes explained 55% of the total variance in the dataset,
only three dimensions of the PC transformed points were considered for
further analysis of niche overlaps and niche width. We used the
NicheRover package on R using a bayesian based probabilistic estimation
of niche space, overlaps, and width
(Swanson et al. 2015).
Directional overlaps were computed by a specified number of Monte Carlo
draws from the parameter list (α = 0.95) that indicated the probability
of finding ‘species A’ in the niche space of ‘species B .’The niche
width was obtainedcalculated by calculating the posterior distribution
of niche size by species (niche size for every posterior sample drawn)
(Swanson et al. 2015).
Niche estimation following the methods described was performed at two
scales (Broad scale: All of Ladakh, three species representing three
niches) and (Regional scales: North-West Ladakh and South-East Ladakh
representing five niches).