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, Reference Schaller and Liu1996). 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, Reference Schaller1998). The species is categorized as Vulnerable on the IUCN Red List (Harris & Leslie, Reference Harris and Leslie2008) 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, Reference Leslie and Schaller2009). Climate change is also expected to affect the long-term availability of suitable habitats for the species (Schaller, Reference Schaller1998), 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., Reference Sharma, Clevers, de Graaf and Chapagain2004; Singh et al., Reference Singh, Yoccoz, Bhatnagar and Fox2009; St-Louis & Côté, Reference St-Louis and 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, Reference Harris and Miller1995; Schaller, Reference Schaller1998); 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, Reference Schaller1998); (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, Reference Jarman1974); and (3) predation risk, here captured by anthropogenic disturbances given the general lack of natural predators of wild yaks in the area (Schaller, Reference Schaller1998), would be a significant factor shaping habitat selection patterns, and wild yaks would avoid areas near human communities (Leslie & Schaller, Reference Leslie and Schaller2009). 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 km2 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, Reference Leslie and Schaller2009). 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, Reference Schaller1998).
Distribution of the wild yak Bos mutus on the Tibetan Plateau, China. The study area covers c. 1.1 million km2, 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., Reference Yu, Xu, Okuto and Luedeling2012) 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, Reference Guisan and Thuiller2005; Austin, Reference Austin2007). The spatial resolution of the environmental variables considered was set to 1 km2. All the candidate variable layers were cropped to the extent of the study area and, if necessary, resampled to a 1 km2 spatial resolution using the nearest neighbour method in the raster library (Hijmans & van Etten, Reference Hijmans and van Etten2012) 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., Reference Hijmans, Cameron, Parra, Jones and Jarvis2005) 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., Reference Hao, Ju, Jiang and Zhu2013).
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, Reference Brown1999; Illius & O'Connor, Reference Illius and O'Connor2000). The topographic ruggedness index, a measurement developed by Riley et al. (Reference Riley, DeGloria and Elliot1999) 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, Reference Schaller1998; Xu et al., Reference Xu, Jiang, Li, Guo, Wu and Cai2006; Leslie & Schaller, Reference Leslie and Schaller2009), predation risk for the species is considered to be shaped primarily by human presence and activity (Leslie & Schaller, Reference Leslie and Schaller2009). 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., Reference Sheehy, Miller and Johnson2006). 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., Reference Nogués-Bravo, Araújo, Errea and Martínez-Rica2007). The melting ice and snowpack provide seasonal fresh water and soil moisture critical to local vegetation communities (Schaller, Reference Schaller1998). 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, Reference Pettorelli2013), 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, Reference Achard and Estreguil1995), 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., Reference Garonna, Fazey, Brown and Pettorelli2009, 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, Reference Guisan and Zimmermann2000). 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, Reference Guisan and Zimmermann2000). 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, Reference Yee and Mitchell1991), maximum entropy (MaxEnt; Elith et al., Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011) and random forests (Breiman, Reference Breiman2001). All models were developed in R using the package biomod2 v. 3.1–48 (Thuiller et al., Reference Thuiller, Georges and Engler2014).
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, Reference O'Brien2007). We excluded some candidate variables to mitigate the effects of inflation caused by the high correlations amongst the predictor variables (Dormann et al., Reference Dormann, Elith, Bacher, Buchmann, Carl and Carré2013).
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., Reference Araújo, Pearson, Thuiller and Erhard2005). 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, Reference Austin2007). When performing MaxEnt, species prevalence was set to 0.1 (Elith et al., Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011). The maximum decision tree of the random forests model was set to 500 (Cutler et al., Reference Cutler, Edwards, Beard, Cutler, Hess, Gibson and Lawler2007). We assessed model performance using three evaluation methods: Kappa (Cohen, Reference Cohen1960), the true skill statistic (TSS; Allouche et al., Reference Allouche, Tsoar and Kadmon2006) and area under the curve (AUC; Swets, Reference Swets1988). Model predictions were categorized as excellent for Kappa > 0.75 (Fleiss, Reference Fleiss1981), TSS > 0.8 (Thuiller et al., Reference Thuiller, Lafourcade, Engler and Araújo2009) or AUC > 0.90 (Swets, Reference Swets1988).
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., Reference Allouche, Tsoar and Kadmon2006; Lobo et al., Reference Lobo, Jiménez-Valverde and Real2008). Variable importance was estimated using a variable permutation algorithm (Breiman, Reference Breiman2001). 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, Reference Hollander and Wolfe1973).
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).
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.
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.
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, Reference Schaller1998) and are therefore not expected to select forage on the basis of quality over quantity (Jarman, Reference Jarman1974). 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, Reference Pettorelli2013, 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., Reference Pettorelli, Ryan, Mueller, Bunnefeld, Jędrzejewska, Lima and Kausrud2011). 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, Reference Harris and Miller1995), 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, Reference Schaller1998). 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., Reference Yu, Xu, Okuto and Luedeling2012). These meadows are associated with high levels of vapour loss (Körner, Reference Körner1999) 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, Reference Schaller1998); 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., Reference Singh, Yoccoz, Lecomte, Côté and Fox2010) and kiang Equus kiang (St-Louis & Côté, Reference St-Louis and 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, Reference Frid and Dill2002; Creel et al., Reference Creel, Winnie, Maxwell, Hamlin and Creel2005). 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, Reference Jarman1974; Clutton-Brock, Reference Clutton-Brock1989). 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, Reference Leslie and Schaller2009).
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., Reference Syfert, Smith and Coomes2013). 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, Reference Pearson and Dawson2003); 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, Reference Pearson and Dawson2003). 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, Reference Kie and Bowyer1999). 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., Reference Araújo, Pearson, Thuiller and Erhard2005). 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, Reference Illius and O'Connor2000). 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, Reference Hardin1960; Leslie & Schaller, Reference Leslie and Schaller2009). Moreover, an increase in hybridization could heighten genetic contamination of wild populations (Leslie & Schaller, Reference Leslie and Schaller2009).
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.
Acknowledgements
Special thanks are due to George Schaller, Madhu Rao, Manuel Mejia, Marcus Rowcliffe, Joel Berger, Clare Duncan, Xueyan Li, Josephine Arthur and two anonymous reviewers for their generous help in developing this study, and valuable comments on the manuscript.
Biographical sketches
Xuchang Liang and Aili Kang are conservation biologists. They have spent years monitoring yaks and studying their ecology on the Tibetan plateau. Nathalie Pettorelli heads the Environmental Monitoring and Conservation Modelling team at the Institute of Zoology, UK, which is interested in developing tools and methodologies to support the sustainable management of natural resources.
