Understanding habitat selection of the Vulnerable wild yak Bos mutus on the Tibetan Plateau

Abstract

We tested a series of hypotheses on drivers of habitat selection by the Vulnerable wild yak Bos mutus, combining distribution-wide sighting data with species distribution modelling approaches. The results indicate that climatic conditions are of paramount importance in shaping the wild yak's distribution on the Tibetan Plateau. Habitat selection patterns were seasonal, with yaks appearing to select areas closer to villages during the vegetation-growing season. Unexpectedly, our index of forage quantity had a limited effect in determining the distribution of the species. Overall, our results suggest that expected changes in climate for this region could have a significant impact on habitat availability for wild yaks, and we call for more attention to be focused on the unique wildlife in this ecosystem.

Introduction

The wild yak Bos mutus is a rare yet iconic large herbivore inhabiting the Tibetan Plateau. One of the largest bovids, wild yaks are also the largest native species in their range, which formerly included China (Gansu, Sichuan, Xinjiang, Tibet, Qinghai), northern India (Ladakh), and Nepal (Schaller & Liu, 1996). The species declined in the 20th century, mainly as a result of excessive hunting; the total number of mature individuals was last estimated to be c. 15,000, in 1995 (Schaller, 1998). The species is categorized as Vulnerable on the IUCN Red List (Harris & Leslie, 2008) and most of the remaining individuals are found in isolated and fragmented populations in the central and northern parts of the Tibetan Plateau. Remnant populations face escalating threats from anthropogenic activities, such as increasing competition with livestock for good grazing areas, and infrastructural development that causes degradation of their habitats (Leslie & Schaller, 2009). Climate change is also expected to affect the long-term availability of suitable habitats for the species (Schaller, 1998), although there is little quantified and spatially explicit information available to inform discussions on potential management options. More broadly, quantitative information on the factors driving patterns in the seasonal distribution of wild yaks is still rare. Studies on Tibetan herbivore species rarely include data from the species’ entire distribution range (Sharma et al., 2004; Singh et al., 2009; St-Louis & Côté, 2014), which prevents the identification of environmental management actions to alleviate further pressures on populations at the necessary scale for the conservation of large species.

We aim to fill this knowledge gap by combining advances in species distribution modelling with sighting data collected across most of the known distribution range of the wild yak. We predicted that (1) the species would exhibit distinct habitat selection patterns between seasons, which has been suggested previously but has not been assessed in a quantitative manner (Harris & Miller, 1995; Schaller, 1998); in particular, we expected preferred habitats of wild yaks during the vegetation growing season to be found at higher altitudes, in more rugged terrain, and closer to glaciers (Schaller, 1998); (2) the species would select for quantity over quality of forage at the distribution-range scale, given that wild yaks are non-selective grazers (Jarman, 1974); and (3) predation risk, here captured by anthropogenic disturbances given the general lack of natural predators of wild yaks in the area (Schaller, 1998), would be a significant factor shaping habitat selection patterns, and wild yaks would avoid areas near human communities (Leslie & Schaller, 2009). The knowledge derived from this study will be used to predict seasonal habitat availability in the context of climate change, which will help to highlight future global conservation challenges on the Tibetan Plateau.

Study area

The study area (Fig. 1) covers c. 1.1 million km² on the Tibetan Plateau. It encompasses the entire Tibet Interior region, from Kunlun in the north to the Gangdise and Nyainqentanglha ranges in the south, with a slight eastward extension to incorporate part of the Sanjiangyuan region in the Qinghai province of China. This region includes most of the known current distribution range of the wild yak (Leslie & Schaller, 2009). Mean annual precipitation follows a decreasing gradient from east to west and from south to north, ranging from c. 500 mm in the south-east to < 50 mm in the north-west. Mean annual temperatures vary from 0 to −6°C, with winter extremes of < −40°C. The Tibetan Steppe is the main ecoregion in the study area. Vegetation is sparsely distributed in alpine meadows, alpine steppes, semi-arid steppes and cold deserts (Schaller, 1998).

Fig. 1 Distribution of the wild yak Bos mutus on the Tibetan Plateau, China. The study area covers c. 1.1 million km², encompassing the entire Tibet Interior region defined by the Kunlun mountains in the north and the Gangdise and Nyainqentanglha ranges in the south, with slight eastward extension to incorporate part of the Sanjiangyuan region in Qinghai Province.

Methods

Presence data

Presence data were collected by the Wildlife Conservation Society (WCS) and partners in 2006, 2008, 2009, 2011, 2012 and 2013. Most of the surveys were conducted within areas known to hold wild yaks; however, the surveys were not primarily designed to collect information on wild yaks, and sightings were opportunistic. Sightings were georeferenced by trained staff, following the field protocol established by WCS China. Ancillary data (e.g. whether sighting data were collected while in vehicle or on foot, number of observers, survey effort) were not systematically collected and therefore could not be taken into account in subsequent analyses. Vehicle surveys were not based on existing road systems but survey effort was shaped by the local topography as well as the distribution of seasonal rivers. While conducting surveys the speed of the vehicle was required to be < 20 km per hour to avoid disturbing wildlife. The total number of independent records of wild yaks was 755; 569 of these sightings were recorded during the non-growing season (October–March; Yu et al., 2012) and the remainder (n = 186) were during the vegetation growing season (April–September).

Environmental variables

We adopted a methodological framework that categorized relevant environmental variables according to limiting factors (i.e. climatic and topographical factors), disturbance (i.e. anthropogenic influence) and resource distribution (i.e. availability of forage and fresh water) (Guisan & Thuiller, 2005; Austin, 2007). The spatial resolution of the environmental variables considered was set to 1 km². All the candidate variable layers were cropped to the extent of the study area and, if necessary, resampled to a 1 km² spatial resolution using the nearest neighbour method in the raster library (Hijmans & van Etten, 2012) in R v. 3.0.2 (R Development Core Team, 2014).

Climate

The 19 Bioclim variables (representative of the years 1950–2000) from the WorldClim dataset (version 1.4; Hijmans et al., 2005) were used to capture current climatic conditions in the study area. To predict future trends in habitat availability we downloaded the Bioclim layers for 2070 under two representative concentration pathways, RCP26 and RCP85, which describe possible climate futures, depending on various trajectories of greenhouse gas concentration. The climate projections used in this study were derived from the HadGEM2-ES climate model, an updated version of the HADGEM model, which has been reported to adequately predict the Tibetan climate (Hao et al., 2013).

Topography

Topography is known to cause variation in the quantity and quality of forage available for large herbivores, and to shape local predation risk (Brown, 1999; Illius & O'Connor, 2000). The topographic ruggedness index, a measurement developed by Riley et al. (1999) to quantify the total altitudinal change across a given area, was calculated based on the downloaded digital elevation model layer GTOPO30 from the U.S. Geological Survey's Long Term Archive (USGS, 2015). Calculations were performed in QGIS v. 2.2.0 (Quantum GIS Development Team, 2014).

Anthropogenic influence

Although natural predators of wild yaks exist on the Tibetan Plateau (e.g. Schaller, 1998; Xu et al., 2006; Leslie & Schaller, 2009), predation risk for the species is considered to be shaped primarily by human presence and activity (Leslie & Schaller, 2009). The distribution of human communities within the study area is relatively dense south of 33°N, and sparse to the north. Livestock rearing is the predominant livelihood. Long-distance nomadism is now rare and pastoral activities normally take place in designated grazing areas near villages (Sheehy et al., 2006). The linear distance between the centre of any given pixel and the nearest village was used as a proxy for anthropogenic disturbance and was calculated for all pixels. Calculations were performed in QGIS using the Proximity function. The shapefile detailing the distribution of villages in the area was provided by WCS China.

Availability of fresh water

Glaciers have important effects on the hydrological cycle of high-altitude regions (Nogués-Bravo et al., 2007). The melting ice and snowpack provide seasonal fresh water and soil moisture critical to local vegetation communities (Schaller, 1998). The linear distance between the centre of any given pixel and the nearest glacier was therefore estimated for all pixels, using the Proximity function in QGIS. The shapefile of glacier distribution was acquired from the GLIMS Glacier Database (GLIMS & National Snow and Ice Data Center, 2005).

Forage

The normalized difference vegetation index (NDVI), one of the most intensely studied and widely used vegetation indices (Pettorelli, 2013), was considered a proxy for forage availability. MODIS Terra NDVI products (MOD13A2, monthly data for 2001–2013) were downloaded using the MODIS Reprojection Tool Web Interface (USGS, 2010). As the reflected light waves captured by satellite sensors can be influenced by a variety of natural phenomena (Achard & Estreguil, 1995), the downloaded data layers were processed in R to convert all negative values to zeros and adjust the anomalous values, which were assumed to reflect atmospheric noise in the MOD13A2 dataset (see Garonna et al., 2009, for full methdology).

Modelling approach

Species distribution models are numerical tools to assist in quantifying species–environment relationships. Increasingly, they are used to gain ecological insights and predict species distributions at large spatial scales (Guisan & Zimmermann, 2000). There are various types of species distribution models that can be used in combination with presence data to assess habitat suitability; the predictive power of a given modelling approach may be context-specific and may vary depending on the study area, variables and resolution considered, as well as the amount of presence data available (Guisan & Zimmermann, 2000). To overcome uncertainties linked to the choice of species distribution model we used three analytical approaches that have been widely employed in species distribution modelling exercises, namely, generalized additive models (Yee & Mitchell, 1991), maximum entropy (MaxEnt; Elith et al., 2011) and random forests (Breiman, 2001). All models were developed in R using the package biomod2 v. 3.1–48 (Thuiller et al., 2014).

Firstly, we explored the importance of the climatic and topographical variables in shaping the current distribution of the wild yak. Because we expected habitat selection to be seasonal, models were run 40 times each for the growing (G_I) and non-growing (NG_I) seasons (Table 1). Secondly, we considered current yak distribution as a function of climatic conditions, topographical factors, availability of forage, glacier distribution and anthropogenic influences. Again, these models were run for the growing (G_II) and non-growing (NG_II) seasons (Table 1). Multicollinearity was checked using the variance inflation factor analysis (R library usdm; O'Brien, 2007). We excluded some candidate variables to mitigate the effects of inflation caused by the high correlations amongst the predictor variables (Dormann et al., 2013).

Table 1 Predictor variables used in this study, with definition, group, and range of values.

*The groups G_I (growing season) and NG_I (non-growing season) included only topographical and climatic variables; G_II and NG_II also included variables capturing information on anthropogenic influence, glacier distribution, and forage availability.

Yak presence data were split independently into 70% for training and 30% for testing (Araújo et al., 2005). Ten thousand background points (representing pseudo-absence for generalized additive models) were selected at random throughout the study area. The generalized additive model was set with four degrees of freedom for smoothing (Austin, 2007). When performing MaxEnt, species prevalence was set to 0.1 (Elith et al., 2011). The maximum decision tree of the random forests model was set to 500 (Cutler et al., 2007). We assessed model performance using three evaluation methods: Kappa (Cohen, 1960), the true skill statistic (TSS; Allouche et al., 2006) and area under the curve (AUC; Swets, 1988). Model predictions were categorized as excellent for Kappa > 0.75 (Fleiss, 1981), TSS > 0.8 (Thuiller et al., 2009) or AUC > 0.90 (Swets, 1988).

Predictions of species presence probability from the best-performing model were converted to presence–absence predictions using a transforming threshold selected as the one that maximized TSS scores (Allouche et al., 2006; Lobo et al., 2008). Variable importance was estimated using a variable permutation algorithm (Breiman, 2001). Information on altitude, terrain ruggedness, and distance to the nearest village and glacier was extracted from all predicted presence pixels for both seasons under the best model of current habitat suitability distribution for wild yaks; these values were then compared between seasons using Wilcoxon one-tailed sum rank tests (Hollander & Wolfe, 1973).

Results

Random forests models generally outperformed general additive and MaxEnt models (Table 2). In accordance with our first prediction, wild yaks showed distinct seasonal patterns of habitat selection; climatic conditions were strong determinants of these patterns at the spatial scale considered (Fig. 2). During the growing season wild yaks appeared to select areas with low levels of fluctuations in monthly precipitation; they also appeared to favour areas with relatively abundant precipitation in the peak summer month (i.e. July). During the non-growing season drier areas with greater fluctuation in monthly precipitation and less extreme winter temperatures were more likely to be preferred (Fig. 2). Preferred habitats during the growing season were found at higher altitudes (W = 2061875023, P < 0.001), closer to glaciers (W = 344529388, P < 0.001) and in more rugged terrain (W = 1800769226, P < 0.001) than those used during the non-growing season (Supplementary Table S1). Contrary to our third hypothesis, however, wild yaks tended to be found closer to villages during the growing season than during the non-growing season (W = 716972327, P < 0.001). All the variables based on NDVI were comparatively of much lower importance in defining habitat selection patterns than the top climatic variables (Fig. 2).

Fig. 2 The importance of individual variables (Table 1) in predicting the distribution of wild yaks under the best model, for (a) the vegetation growing season and (b) the non-growing season. The best model was run 40 times for each season, and variable importance was evaluated independently for each run.

Table 2 Performance of random forests, general additive, and MaxEnt models in terms of area under the curve (AUC), true skill statistic (TSS) and Kappa, for the growing and non-growing seasons. Each model was run 40 times for each season, and model performance was evaluated independently for each run.

Based on these results it is likely that under the RCP26 scenario the distribution of suitable habitats for wild yaks would expand by 146 and 35% by 2070 in the growing and non-growing seasons, respectively. Under the RCP85 scenario, however, the distribution of suitable habitats during the growing season would expand by 194%, whereas the availability of suitable habitats during the non-growing season is expected to decrease by 76% (Fig. 3). Shifts in the distribution of suitable habitats are also expected to occur. Based on our analyses the distribution of suitable habitats during the growing season could shrink by 69% under RCP26, and 74% under RCP85. Likewise, the distribution of suitable habitats during the non-growing season could shrink by 49% under RCP26, and 98% under RCP85 (Supplementary Table S2). For information on topographical features of potential suitable habitats under both RCPs, see Supplementary Table S3.

Fig. 3 Distribution of suitable habitat for wild yaks in (a) the vegetation growing season and (b) the non-growing season, and the predicted distributions under the RCP26 scenario (a1 & b1) and under the RCP85 scenario for both seasons (a2 & b2).

Discussion

Our results largely support current expectations about the factors shaping wild yak distribution on the Tibetan Plateau, showing that the species’ habitat selection patterns are seasonally distinct and largely driven by climatic factors. Yet two of our predictions were not well supported by our findings. The first pertains to the importance of forage quantity in driving habitat selection. Wild yaks are non-selective grazers (Schaller, 1998) and are therefore not expected to select forage on the basis of quality over quantity (Jarman, 1974). Although we expected forage biomass to be a key factor in determining wild yak occurrence, our results show that most variables based on NDVI play little if any role in shaping wild yak distribution. Unlike previously reported cases where NDVI could be linked to the distribution of large herbivores (see Pettorelli, 2013, for a review), variables based on this index may not have captured vegetation dynamics correctly in our study area as a result of high soil reflectance (Pettorelli et al., 2011). However, our results may also suggest that wild yaks select for quality over quantity of forage to an extent beyond our initial expectation. The nutritious Kobresia-dominant moist meadows, favoured by wild yaks in summer according to empirical observations (Harris & Miller, 1995), are not as productive in terms of vegetation biomass as other vegetation types, such as Stipa grasslands, which are more widely distributed in our study area (Schaller, 1998). Pixels with higher NDVI values would thus fail to capture the distribution of these favoured, yet less productive, meadows. The low level of fluctuation in monthly precipitation, and the abundant precipitation in July (the two conditions identified as being key to capturing wild yak distribution during the growing period) are also key factors determining the biomass and nutrient value of Kobresia-dominant moist meadows (Yu et al., 2012). These meadows are associated with high levels of vapour loss (Körner, 1999) and are therefore dependent on water availability to prevent desiccation. In July, in particular, vegetation in the Kobresia-dominant moist meadows is normally at its early phenological stages (Schaller, 1998); timely and abundant precipitation could thus be particularly beneficial to plant development in these meadows. Studies from other parts of the plateau on the Tibetan argali Ovis ammon hodgsoni (Singh et al., 2010) and kiang Equus kiang (St-Louis & Côté, 2014) similarly suggest that forage quality can be a key factor shaping habitat selection patterns for these large herbivores. At this stage it is difficult to be conclusive about the role of forage quantity and quality in driving wild yak habitat selection; further research is needed.

The second prediction that our results failed to support is that wild yaks avoid human settlements, especially during the period when forage is relatively abundant and when there is thus no need to take greater risks associated with proximity to humans (Frid & Dill, 2002; Creel et al., 2005). The low influence of anthropogenic disturbance on the distribution of wild yaks may suggest that individuals in the area are essentially unaffected by human distribution during the growing season, but this result may also be underpinned by the spatial proximity of villages to Kobresia-dominant moist meadows. Another potential explanation is based on the distribution of domestic yaks, found near villages. Habitat selection patterns of polygynous male herbivores are likely to be dependent on the spatio–temporal distribution of females during the mating season (Jarman, 1974; Clutton-Brock, 1989). One could expect wild males to be attracted by the presence of a large number of domestic females, with no competitors present. This hypothesis is supported by the increasingly reported mingling and hybridization of wild and domestic yaks in Tibet (Leslie & Schaller, 2009).

There are a number of caveats associated with our data and modelling work. Firstly, apart from geographical coordinates and group size the survey teams recorded no other observations (e.g. topography, climatic conditions, primary productivity). All the environmental information used in the analyses was therefore derived from global products, which have not been validated locally. Future research should ground-truth these products to ascertain the robustness of our conclusions. Secondly, our proxy of anthropogenic disturbance does not differentiate disturbance resulting from human presence from disturbance resulting from livestock. This lack of differentiation is a result of a lack of information on the spatial distribution of people and livestock in the area. As efforts to gather such data increase it would be interesting to contrast the influence of humans and livestock on the distribution of wild yaks. Thirdly, the dataset may have been biased by the survey methods. In the growing season, in particular, limited accessibility to various areas can limit survey efforts to regions closer to villages, and therefore our dataset may not capture the full range of environmental conditions where yaks can be found during this period. This sampling bias could lead to underestimation of the distribution and size of suitable habitats during the growing season, as well as misidentification of the ecological forces shaping the distribution of the species (Syfert et al., 2013). Moreover, this study, based on a series of correlative modelling approaches, intrinsically assumes that wild yaks are living in equilibrium with their environment (Pearson & Dawson, 2003); however, the observed distribution may not reflect the optimal patterns of habitat selection but rather habitat use as constrained by a number of factors, including those associated with the presence of livestock. To address this, future large-scale studies should attempt to incorporate information on the distribution of domestic yaks when modelling the distribution of wild yaks. Various biotic interactions, as well as the dispersal ability of wild yaks, need to be taken into account to identify scale-dependent limiting factors and consequent patterns in habitat selection (Pearson & Dawson, 2003). Another limitation to this study is that we did not consider the influence of sex. Dimorphic ruminants can be substantially divergent in their niche requirements (Kie & Bowyer, 1999). We were unable to explore differences in habitat selection patterns between males and females because the sex of the individuals sighted was not recorded reliably. The seasonal patterns identified should thus be understood as averaged results based on the unknown sex ratio. Finally, uncertainties associated with the modelling approaches considered should be acknowledged (Araújo et al., 2005). Predictions derived from these models vary substantially; for example, if we adopted the predictions of the generalized additive model (which is also acceptable in terms of AUC and TSS) the importance of factors such as altitude and mean temperature in determining suitable habitats for wild yaks in summer would be much higher than suggested by the random forests model (Supplementary Fig. S1), and suitable habitats would be more extensive in both seasons (Supplementary Table S2). Such method-induced variability highlights the importance of interpreting model outputs with caution.

An important contribution of this study is its quantification of the possible impact of climate change on the availability of suitable habitats for wild yaks. According to our current knowledge wild yaks are mostly found at 33–36°N, in regions that are likely to be severely affected by climate change. In terms of conservation priorities for the species, its autumn and winter habitats appear to be more vulnerable to the impacts of climate change than its spring and summer habitats, and a lower winter to summer habitat ratio may represent a high risk to population stability. In a seasonal grazing system, ungulate population size can be more sensitive to limited availability of resources during the non-growing season than to abundance of resources during the growing season (Illius & O'Connor, 2000). In the future the distribution of suitable habitat during the growing season is more likely to be threatened by anthropogenic activities than by climate change. Any increase in the distribution of suitable habitats will create economic opportunities for domestic yak herders, which could result in resource competition between wild and domestic yaks at local scales, and increase the potential for disease transmission between groups (Hardin, 1960; Leslie & Schaller, 2009). Moreover, an increase in hybridization could heighten genetic contamination of wild populations (Leslie & Schaller, 2009).

Our results suggest that increasing dispersal opportunities for local yak populations should be a key component of any conservation scheme aiming to mitigate the impact of climatic change, helping them to track shifting climatic zones and colonize new territories. Our findings also suggest that the number of domestic yak holdings should be more strictly controlled in communities adjacent to known wild yak populations. Livestock grazing activities should be limited to designated areas to minimize competition for the winter resources of wild yaks. These two points are especially relevant for two regions that include parts of the Ali and Naqu prefectures, where there are high densities of wild yaks. These regions are likely to remain suitable for the species under both RCP scenarios during the growing season; however, they may not remain so during the non-growing season. They are not within the borders of any protected area and are subject to high levels of human activity. Conservation interventions may be necessary in these regions and we recommend the establishment of monitoring systems as soon as possible, to assess any direct threats, such as illegal hunting. Current patterns of land use (e.g. grazing sites for domestic yaks) should be evaluated and possibly rearranged to take into consideration the habitat needs of wild yaks. Lastly, we recommend the rapid definition and implementation of a plan to connect these regions to the nearest protected areas that contain other wild yak populations.