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:
  1. 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.
  2. 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.
  3. 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.
  4. 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).