A diet rich in C3 plants reveals the sensitivity of an alpine mammal to climate change

Abstract Plant–herbivore interactions provide critical insights into the mechanisms that govern the spatiotemporal distributions of organisms. These interactions are crucial to understanding the impacts of climate change, which are likely to have an effect on the population dynamics of alpine herbivores. The Royle's pika (Ochotona roylei, hereafter pika) is a lagomorph found in the western Himalaya and is dependent on alpine plants that are at risk from climate change. As the main prey of many carnivores in the region, the pika plays a crucial role in trophic interactions. We examined topographical features, plant genera presence and seasonal dynamics as drivers of the plant richness in the pika's diet across an elevational gradient (2,600–4,450 m). We identified 79 plant genera in the faecal pellets of pikas, of which 89% were forbs, >60% were endemic to the Himalaya, and 97.5% of the diet plant genera identified followed the C3 photosynthetic pathway. We found that, during the premonsoon season, the number of genera in the pika's diet decreased with increasing elevation. We demonstrate that a large area of talus supports greater plant diversity and, not surprisingly, results in higher species richness in the pika's diet. However, in talus habitat with deep crevices, pikas consumed fewer plant genera suggesting they may be foraging suboptimally due to predation risk. The continued increase in global temperature is expected to have an effect on the distribution dynamics of C3 plants and consequently influence pika diet and distribution, resulting in a significant negative cascading effect on the Himalayan ecosystem.

Determining the diet of a herbivore is fundamental to understanding trophic interactions and assessing dietary plasticity to climate change, and for developing effective monitoring, management and conservation strategies (Bernstein et al., 2007).
DNA metabarcoding using high-throughput sequencing enables the identification of dietary species using DNA extracted from faecal samples more accurately than traditional visual and microscopic observation of faecal samples (Kartzinel et al., 2015). This technique has been used to identify herbivore gut contents, revealing cryptic functional diversity and niche partitioning (Kress, García-Robledo, Uriarte, & Erickson, 2015). Plastid genes such as rbcL (Group et al., 2009) and nuclear ribosomal internal transcribed spacer (ITS; Hollingsworth, 2011) are commonly used for plant metabarcoding (Hollingsworth, 2011). We aimed to analyse faecal pellets using metabarcoding to identify the plants in the diet of the Royle's pika and examine the effects of talus characteristics, topography and plant richness; abundance; and seasonal dynamics (pre-and postmonsoon) on diet, across five sites in the western Himalaya, India. Food availability (plant species' presence and abundance ;Huntly, Smith, & Ivins, 1986;Dearing, 1995Dearing, , 1996Wilkening et al., 2011;Bhattacharyya et al., 2013Bhattacharyya et al., , 2014aBhattacharyya & Ray, 2015), habitat topography (elevation, slope, aspect;Walker, Halfpenny, Walker, & Wessman, 1993;Deems, Birkeland, & Hansen, 2002;Wilkening et al., 2011;Rodhouse et al., 2010;Gurung et al., 2017) and predation risk (rock cover, crevice depth, distance to the nearest area of talus; Calkins, Beever, Boykin, Frey, & Andersen, 2012;Bhattacharyya et al., 2014a,b;Castillo, Epps, Davis, & Cushman, 2014; significantly impact the foraging ecology of the Royle's pika and potentially influence access to nutritive plants, and thereby affect individual fitness (Bhattacharyya, 2013). For talus-dwelling Ochotona spp., talus size (area) and connectivity between talus are known to influence the habitat occupancy (Franken & Hik, 2004). In Royle's pikas, a high proportion of rock cover in talus habitat provides refuge from predation risk, which in turn increases their habitat occupancy and abundance (Bhattacharyya et al., 2014a. Topographical features determine the distribution and abundance of alpine plants (Bruun et al., 2006). The prehistoric distribution of small mammals (Ochotonidae), highly adapted to arctic or alpine environment, was closely associated with preferred C 3 food plant diet which is high in protein and moisture content (Ge et al., 2012). Whilst the interplay between plants and herbivores is a key determinant of community structure (see Dolezal et al., 2016;Wisz et al., 2013), climate-induced expansion of C 4 plants is believed to have resulted in extinction and range contraction in pikas (Ge et al., 2012). Furthermore, in highly seasonal montane habitats, dietary constraints in herbivores tend to be strongly linked to quality of forage available. In the Himalayan region, both plant quality and available biomass may act as constraints for pikas. Seasonal dynamics in diet selection can reflect dietary adaptations (plasticity) in a seasonal alpine habitat. Our Royle's pika diet analysis will provide insights into the plant genera selected during foraging and test the following hypotheses: 1. Seasonal difference in diet (premonsoon versus postmonsoon) will demonstrate dietary flexibility, and the pikas will select plants with the highest nutrient gain possible from the feeding habitats available.
2. Larger areas of talus habitat will result in more plant availability and will increase species richness of the diet.
3. Talus characteristics such as rock cover and depth of crevices would increase species richness in the diet by providing refuge BHATTACHARYYA ET AL. | 251 from predators and allowing pikas to access larger foraging ground.
4. Forbs and grasses constitute the largest proportion of their diet. 5. Pikas prefer diets rich in C 3 plants, which signify potential threat due to climate change as C 3 plants are dependent on high rainfall and low temperature.

| Study area, collection of faecal pellets and habitat data
The Royle's pika inhabits rock crevices in talus fields with a home range of approximately 50 m 2 Kawamichi, 1968  Following , we selected 50 m 2 survey plots in each talus habitat, and these were searched for fresh piles of faecal pellets (moist, dark brown/black). One site (

| DNA metabarcoding of faecal pellets
The Royle's pika's diet is rich in secondary metabolites , which often inhibit downstream enzymatic reactions in polymerase chain reaction (PCR; Weishing, Nybom, Wolff, & Meyer, 1995). Therefore, we homogenized 20-30 mg of the faecal samples and extracted DNA using QiAamp DNA Stool Kit (Qiagen, Manchester, UK) following minor modifications in the manufacturer's protocol (e.g., overnight incubation at 56°C with ASL stool lysis buffer; Qiagen Inc., Germany). As Royle's pika were the only lagomorph species present in the study area (Bhattacharyya, 2013;Green, 1985) with very distinct faecal pellets, the chances of misidentification of the faecal samples were ruled out.
We amplified the ITS2 region of plant nuclear DNA using primer  Levin et al., 2003;Kress et al., 2009;Yoccoz et al., 2012). Primers had overhang adapter sequences added to the 5′ end for the initial PCR amplification, following Campbell, Harmon, and Narum (2015

| Identification of plant diet from sequencing data
The bioinformatic analyses were performed using Iceberg, the High Performance Computing Cluster at the University of Sheffield, UK. The paired-end reads were filtered for quality (minimum quality score 20 over a 4 bp sliding window) and any Illumina adapter sequences removed using Trimmomatic v 0.32 (Bolger, Lohse, & Usadel, 2014), retaining only reads of at least 90 bp in length. For ITS2, filtered sequences were then aligned using FLASH (Magoč & Salzberg, 2011), aligned sequences with matches to the ITS2 primer sequences only were extracted, and primer sequences were removed using the "trim_seqs" command in Mothur BHATTACHARYYA ET AL.
| 253 (Schloss et al., 2009 The "derep_fulllength" and "uchime2_denovo" commands were used in USEARCH software v 9.2.64 (Edgar, 2010) to eliminate all sequences which had less than 10 copies per sample and any chimeric sequences. Each unique ITS2 sequence found in the data set was then compared against the NCBI GenBank nucleotide database using the BLAST algorithm (Altschul et al., 1997) to assign a taxonomic unit for each plant sequenced identified. Only matches with at least a 97% identity to the reference sequence were retained for downstream analysis. The software METAGEN-OME ANALYZER v 4 (MEGAN, Huson et al., 2016) was used against NCBI taxonomic framework for mapping and visualization of the BLAST results, keeping all default parameters for the LCA assignment algorithm except the bit score minimum support threshold, which was set at 1%.
For rbcL, the analysis pipeline differed slightly as we did not expect a

| Information of ecological and evolutionary linkages of plants in Royle's pika diet
The adaptation ability, ecological requirements, physiology and nutri-

| Plant composition and seasonal dynamics in diet
We compared the plant diet composition at family and genera level. We used multiple response permutation procedure (MRPP; Mielke, Berry, & Johnson, 1976;Zimmerman, Goetz, & Mielke, 1985) to understand seasonal variation in diet composition at plant genera level across multiple sites. We used contingency table analysis to test for heterogeneity in plant family in diet across sites, which was assessed by G tests followed by partitioned analyses (Sokal & Rohlf, 1995). In addition, we have used generalized linear models (GLMs) to understand the effect of elevation on plant diet richness (number of genera) across seasons and comparison of plant families across sites sampled in the same season. The significance of fixed effects was evaluated with Wald's chi-square tests (Bolker et al., 2009). To understand whether relative contribution of vegetation type (forbs, grass, shrub and tree) in pika diet is proportional to their availability in the environment, we conducted compositional analysis (Aebischer, Robertson, & Kenward,1993) using Adehabitat v 1.8.20 package in R (Calenge, 2006). We restricted our analysis to the vegetation types which qualified within the criteria of a minimum of two data points greater than zero per vegetation category in the environment data set (Calenge, 2006). In addition, the contribution of each plant genus to overall pika diet with their corresponding evolutionary origin and distribution range was visually explored using box plots.

| Modelling effect of talus habitat on diet richness
We investigated the effect of food availability (tree, shrub, grass and forbs cover), predation risk (talus area, distance between talus, and depth of crevices) and talus topography (elevation, aspect and slope) on plant diet richness. We fitted 32 GLM models and used AIC-based multimodel inference to identify well-supported statistical models that describe the relationships between plant diet richness and biological parameters relevant to foraging ecology of Royle's pika (Table S1) including models with ΔAICc < 2 (Burnham & Anderson, 2002). The relative importance (RI) of each parameter after model averaging was calculated by summing Akaike's weight (wi) across all models in which the parameter was present. All data were checked for normality and corrected for overdispersion if required. We tested for possible collinearity of the explanatory variables using Pearson's correlation analysis in the global model; the mean correlation was 0.14, and the strongest was 0.64. Correlated variables were used only in interactive models but not in additive models (Graham, 2003). All analyses were conducted in R v 3.3.3 (The R Foundation for Statistical Computing, http://www.r-projec t.org/).

| Plant composition and seasonal dynamics in pika diet
We successfully obtained information from 110 faecal pellets (after controlling for data quality) representing 66 out of 104 sampled rocky talus plots. A total of 79 plant genera (ITS2 = 54, rbcL =13, 12 were common to both) were identified. We retrieved 62 genera and 28 families of forbs, ten genera and one family of grass, three genera and three families of shrubs, and three genera and three families of trees (Table S2).
Of the 32 plant families recorded in the pikas' diet, Asteraceae, Poaceae, Primulaceae, Ranunculaceae, Rosaceae and Scrophulariaceae showed significant differences across sites in the mean proportion of occurrence within diet (Table 1). We found no significant difference in proportion of plant families in the premonsoon (GLM: In the premonsoon season, we found significant differences between sites in the proportion of plant family per faecal sample (GLM: F 2 = 126.92, p < 0.0001; proportions listed in Table 1)  Forbs constituted 89% of richness, which was significantly higher than both shrubs (7%) and grasses (3%) (GLM: F 2 = 48.91, p < 0.03).

| Effect of talus habitat on plant diet richness
Model-averaged estimates derived from the 90% model set agreed with the best approximating model with two variables: talus area cover and crevice depth, which were detected as significant predictors for diet richness; each having RI values of 1.0. Talus area has significant positive influence on plant diet richness. However, crevice depth showed a negative influence on plant richness (Table 2).

| DISCUSSION
This is the first study of the diet for any Himalayan pika species that uses noninvasive sampling and metabarcoding. It allowed us to quantify genus richness and revealed cryptic aspects of functional diversity and useful insights into niche partitioning across an elevational gradient. Our study provides an excellent example of how DNA metabarcoding can be used in understanding diet and feeding preferences of an elusive herbivore species found in fragmented alpine terrain and is applicable for other herbivores dependent on C 3 and C 4 plants. We found DNA metabarcoding outperformed traditional methods by revealing the huge diversity of plant genera consumed by pikas and their reliance on species endemic to the Himalayas.
Earlier studies based on traditional visual and microscopic BHATTACHARYYA ET AL.

Royle's pika diet
Our study revealed that pikas exhibit dietary flexibility with high genus richness in their premonsoon diet at lower elevation, possibly due to a longer growing period in this habitat. Thick snow cover delays the T A B L E 1 Presence of plant families in the Royle's pika's diet across each sampled site (TUN: Tungnath; MAD: Madmaheshwar; NAN: Bedni-Roopkund; RUD: Rudranath; and HAK: Har ki Doon), with results from G test (G statistics and p values), to estimate differences in the mean proportion of plant families within diet. Proportion of plant family is presented for each site; those highlighted in bold differ across sites at p < 0.05 Note. a Plant family not considered for comparison across sites as they were found only in one study site (in grey shade).

Family Total genus Premonsoon Postmonsoon G p TUN (n = 27) MAD (n = 13) NAN (n = 32) TUN (n = 26) RUD (n = 4) HAK (n = 8)
beginning of the growing period of alpine plants in higher elevation areas (Inouye, 2008). Therefore, plants experience comparatively longer growing periods at lower elevations (Körner, 2003). However, the lack of precipitation and low temperature during the postmonsoon period leads to the end of the growing period for alpine plants. This may explain why no variation in diet plant richness with elevation was  Khlebnikova, 1976;Revin & Boeskorov, 1990), Royle's pikas show a preference for plants such as forbs (e.g., Potentila spp., Primula spp. and Anaphalis spp.).
The contribution of forbs in the pika's diet was found to be proportionately higher than its availability in the environment possibly due to the high nutrient value of forbs compared to other vegetation types, such as grasses, shrubs and trees .
Talus-dwelling Ochotona species usually exhibit prominent hay building activity in the summer to cache food for survival during the winter (Bhattacharyya & Ray, 2015;Hudson, Morrison, & Hik, 2008).
This food caching helps in lowering the quantity of secondary metabolites during winter consumption (Dearing, 1997  Predation risk has significant impacts on prey populations, either by direct predation-mediated mortality or indirectly by altering their physiology and behaviour (Lima, 1998;Lima & Dill, 1990;Schmitz, Beckerman, & O'Brien, 1997;Sinclair & Arcese, 1995). Large talus areas with high rock cover and vegetation help in forage selection by reducing other constraints such as energy demands (Stephens & Krebs, 1986).
Large rock talus fields with crevices possibly act as escape cover from predators and allow pikas to access extended foraging grounds and diverse array of food plants, and hence have a positive influence on diet richness. Furthermore, Royle's pikas are sensitive to high temperature, and stable microclimate talus habitat serves as a refuge from harsh climate as well as predators (Bhattacharyya et al., 2014b). Royle's pika utilizes crevices in their talus habitat to build their nests . The structure of rock talus and availability of small crevices (< 15 cm) appear to govern occupancy of Royle's pika habitat as it reduces predation risk; wide crevices probably make it easier for predators such as weasels and red foxes to catch pikas . We found deep crevices had a negative effect on plant diet richness. The negative effect of deep crevices probably indicates foraging under fear response. The deep crevices often also have wider openings that increase predation risk from small size predators such as weasels and red foxes .  Holmes, 1991;Morrison et al., 2004;Roach, Huntly, & Inouye, 2001).

Royle's pika diet
The majority (97.5%) of diet plants in this study were C 3 food plants, The bold value indicates statistically significant p value which is less than 0.05.
T A B L E 2 Model-averaged predictive models for effect of each parameter on plant diet richness with AIC < 2 and relative importance (RI) of each parameter has shown that expansion of C 4 plants with low nutrient availability might have led to the extinction and contraction of distributional range of herbivorous mammals which were dependent on C 3 plants during late Miocene (Cerling, Ehleringer, & Harris, 1998;Ehleringer, Cerling, & Dearing, 2002;Ge et al., 2012;MacFadden & Ceding, 1994;Osborne, 2008;Osborne & Beerling, 2006). The prehistoric distribution of pikas was closely associated with distribution of their preferred C 3 food plants, such as those with high protein and moisture content (Ge et al., 2012). Climate change has also been found to alter the distribution and abundance of alpine plants based on their temperature sensitivity (Scherrer & Körner, 2011). Replacement of C 3 plants by C 4 plants across large landscapes during the late Miocene resulted in the extinction of a high number of pika species (Ge et al., 2012). Gottfried et al. (2012) suggested that the thermophylization process could lead to the replacement of cold and moist environment plant species (e.g., C 3 plants) with an abundance of warm and dry environment plant species, (e.g., C 4 and CAM plants).
Given that the Royle's pika diet consists of C 3 plant species (e.g., Potentila, Anaphalis), the increase in temperature and rainfall patterns is bound to influence the distribution range of plant species as well as pikas. In the past 25 years , the Himalayan arc has experienced significant decreases in winter precipitation (17 mm Alpine species adapted to cold climatic conditions are more vulnerable to global warming (Hughes, 2000). Moritz et al. (2008) indicated significant changes in distributions and range contraction in mountain-dwelling small mammals. Walther et al. (2002) predicted upslope range shifts in animals to cope with changing (warming) environmental conditions. However, habitat fragmentation might not allow small mammals such as pikas to move their ranges fast enough to track shifts in suitable microclimates (Ray et al., 2012;Schloss, Nuñez, & Lawler, 2012). Given that the Royle's pika can cope with physiological stress (e.g., hypoxia) in high-elevation environments, the upper limit of the distribution of diet plants (> 4,500 m) could still pose a threat to their survival and fitness. Isolated high altitude mountain habitats, also known as "Sky-islands," differ significantly in environmental conditions from intervening valleys (Shepard & Burbrink, 2008). These intervening valleys were found to have significant impacts on the dispersal of cold-adapted mountain-dwelling species, such as pikas, across sky-islands (Galbreath, Hafner, Zamudio, & Agnew, 2010). Hence, the Royle's pika would potentially need to travel across the lowland valleys to find favourable habitat to cope with the changing environment. Apart from high sensitivity to warm environments (Bhattacharyya et al., 2014b) and predation risk , absence of diet plants below 1,500 m elevation might hamper such contemporary migration through lowland valleys. Furthermore, loss of distribution range and connectivity between habitat patches along with nutritional stress might also make small pika populations more vulnerable to both climatic changes as well as other threats, such as the susceptibility to infectious diseases (Biebach & Keller, 2009;Epstein, 2001;Hanski & Gilpin, 1991;Harvell et al., 2002).
India holds 28% of the flora which is endemic to the Himalayan region (Chitale, Behera, & Roy, 2014). We found a high dependency of the Royle's pika on endemic plants, and 53%-64% of dietary plant genera (e.g., Geum and Fragaria) are endemic to the Himalayan region. Increases in summer temperature and precipitation and the frequency of freeze-thaw cycles are predicted to have detrimental, multidimensional and spatially variable impact on these endemic alpine plants (Dolezal et al., 2016). Recent research has suggested around 23.9% to 41.34% reduction in the distribution of endemic plants in various biodiversity hotspots in Himalaya by 2080 (Chitale et al., 2014). Therefore, changes in the distribution and abundance of Himalayan endemic plants might have significant negative impacts on the overall nutritional ecology of the Royle's pika. Lack of a plant genetic reference library from the Himalayan region, especially from our study site, did not allow us to achieve species-level identification of pika diet. Future research is needed to build a plant genetic resource reference library for insights to be gained at the plant species level.

| CONCLUSION S
This is the first study of the Royle's pika diet using DNA metabarcoding and successfully revealed feeding preferences in this climatesensitive herbivore. We were able to quantify the genus richness in the pika's diet in various habitats and at different elevations. We revealed the high dependency of the pika on C 3 plants and plants endemic to the Himalayan region. We showed that seasonal differences in the diet are associated with elevation and that diet varies at different sites due to differences in topography and climate. A significant amount (> 89%) of the pika's diet consisted of forbs. The genus richness in the diet is strongly predicted by the size of the talus available to the pika for foraging, with large areas resulting in the highest levels of diversity observed in the diet. Crevice depth negatively influences plant diet richness due to increased predation risk. The continued increase in temperature could have a significant effect on the distribution of C 3 plants and reduce the amount of plant species available on which the pika can feed, which would be likely to impact both pika numbers and their distribution, especially in lower elevation habitats.

ACKNOWLEDG EMENTS
This study was partially funded by the Department of Biotechnol- Gayatri Anand for proofreading.

DATA ACCESSIBI LITY
GPS locations of all sample collection points, R code and data for generalized linear modelling, DNA sequences generated during the study can be accessed in the Dryad database (https://doi.org/10. 5061/dryad.nt40m53).

AUTHORS' CONTRIBUTI ON
S.B. collected the samples and performed the laboratory work.