Projected distribution and climate refugia of endangered Kashmir musk deer Moschus cupreus in greater Himalaya, South Asia

Kashmir musk deer Moschus cupreus (KMD) are the least studied species of musk deer. We compiled genetically validated occurrence records of KMD to construct species distribution models using Maximum Entropy. We show that the distribution of KMD is limited between central Nepal on the east and north-east Afghanistan on the west and is primarily determined by precipitation of driest quarter, annual mean temperature, water vapor, and precipitation during the coldest quarter. Precipitation being the most influential determinant of distribution suggests the importance of pre-monsoon moisture for growth of the dominant vegetation, Himalayan birch Betula utilis and Himalayan fir Abies spectabilis, in KMD’s preferred forests. All four Representative Concentration Pathway Scenarios result an expansion of suitable habitat in Uttarakhand, India, west Nepal and their associated areas in China in 2050s and 2070s but a dramatic loss of suitable habitat elsewhere (Kashmir region and Pakistan-Afghanistan border). About 1/4th of the current habitat will remain as climate refugia in future. Since the existing network of protected areas will only include a tiny fraction (4%) of the climatic refugia of KMD, the fate of the species will be determined by the interplay of more urgent short-term forces of poaching and habitat degradation and long-term forces of climate change.

. Occurrences (n = 136) of Kashmir musk deer against the backdrop of precipitation of the driest quarter, which was identified by this study as the most important predictor of musk deer distribution. A total of 123 points were obtained from 10 protected areas situated between central Nepal and Afghanistan through primary data collection efforts between 2003 and 2016. These data were supplemented with 13 points from Afghanistan and Kashmir region obtained from published sources. The map was plotted using R 3.4.3 (R Foundation for Statistical Computing, Vienna, Austria, https:// www.R-project.org).
In this study, we modelled the distribution of the least studied endangered Kashmir musk deer with the following objectives: (1) understand the current and expected future distribution of the species under different climatic scenarios, (2) identify potential climatic refugia for the species, and (3) determine the adequacy of the existing network of protected areas for protecting this species in future.

Results
important predictors of species distribution. Out of 22 selected covariates, we selected for species distribution modeling, four were found to have disproportionately high influence in determining the distribution of the KMD (Fig. 2). These include precipitation of the driest quarter "bio17" (relative influence = 65%) and annual mean temperature "bio1" (relative influence = 18%), as the two top most predictors. Water vapor "vapr" and precipitation during coldest quarter "bio19" played a minor but still important role in species distribution. Collectively, the predictors possessed relative influence equal to 93% in model performance.
projected distribution under current climatic conditions. The prediction of the species distribution models yielded probability of occurrence or habitat suitability (Fig. 3). KMD is predicted to inhabit a belt of high Himalaya that stretches from central Nepal to north-west of India, reaching Afghanistan through the Kashmir Region of India and Pakistan. However, the suitable habitats do not occur in a continuum throughout the high Himalaya. Suitable habitats are located west of Annapurna Himalaya range including Mustang (hereafter Annapurna region), west of Annapurna region (hereafter, west Nepal), the northern part of Uttarakhand state of India (hereafter, Uttarakhand), west of Uttarakhand in Himachal Pradesh (hereafter Himachal), the Kashmir region of India and Pakistan (hereafter, Kashmir region), and along and close to the northeast border between Pakistan and Afghanistan (hereafter, Pak-Afghan border). Additionally, small patch of suitable habitat observed in Lubu area in China nearby north west border between Nepal and China. Few patches of suitable habitat occur in the the eastern half of Nepal and their adjacent area in Tibet as this area supports the closely related Himalayan musk deer.
When the continuous probability surface was classified into categories, we observed that highly suitable (probability >0.7) and suitable (probability 0.5-0.7) habitats are distributed in patches throughout the Himalaya from central Nepal to north east Afghanistan ( Supplementary Fig. S1). The continuous patches of high quality habitat Figure 2. The relative influence of the predictors in Kashmir musk deer distribution; predictors are listed on y-axis. We started our model construction process with 26 potential predictors (Table 1). However, four of these covariates (land cover, aspect, altitude and distance to the nearest water source) were not available for future climate scenarios at our desired spatial resolution or required enormous computational resources. These covariates were some of the least important predictors (relative influence <3), and therefore were dropped from further analysis. The remaining 22 predictors (Table 1)   will differ in four scenarios. In 2050s (RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5), most of the suitable habitat will be available in west Nepal and Uttarakhand as well as their adjacent area in Tibet along the border between three countries (Nepal, India and China). In RCP 2.6 scenario in 2050, some patches of suitable habitat will appear in Himachal between Kashmir region and Uttarakhand in future whereas most of the habitat will disappear from Kashmir region and Pak-Afghan border. In RCP 4.5, RCP 6.0 and RCP 8.5 scenarios, by 2050, KMD's suitable habitat will disappear from Kashmir region, Pak-Afghan border and Annapurna region. In all four scenarios during 2070s, the suitable habitat from Pak-Afghan border, Kashmir region and Himachal region will disappear. In the rest of Himalaya to the west of Uttarakhand (the Indian Himalaya, all Himalaya of Pakistan and Afghanistan Himalaya), currently observed suitable and highly suitable habitats will be entirely lost in future (2050s and 2070s) in three scenarios (RCP 45, RCP 60 and RCP 85). In the Pak-Afghan border, Pakistan, Kashmir region, and central Nepal regions only marginally suitable are predicted to be present in future during 2050s and 2070s (Fig. 4). Compared to the other RCP scenarios (RCP 2.6, RCP 6.0 and RCP 8.5), by 2050 the habitat for KMD under the RCP 4.5 will expand more and the habitat of this deer is expected to be distributed densely within west Nepal, Uttarakhand and their associated areas in Tibet (Fig. 4). By 2070, habitat in RCP 4.5 and RCP 8.5 scenarios will depict similar pattern (Fig. 4). However, the habitat of KMD will be expanded in Uttarakhand, west Nepal and their associated areas in Tibet in RCP 8.5 scenarios (Fig. 4). The habitat in RCP 2.6 in 2070s will be found sparsely in east of Uttarakhand, west of west Nepal and their adjacent area in Tibet compare to other three scenarios. Habitat of KMD in three scenarios (RCP 26, RCP4.5 and RCP 6.0) will decrease in 2070s compare to 2050s in Uttarakhand, west Nepal and their associated areas in Tibet. However, habitat of KMD in RCP 8.5 will increase in 2070s compare to 2050s (Fig. 4).
In future, HSH (Highly Suitable Habitat) for KMD will dominate over SH (Suitable Habitat) and MSH (Marginally Suitable Habitat) in Uttarakhand and west Nepal in all four climatic scenarios. In Tibet, HSH dominate MSH and SH in three cases i.e. RCP 4.5 and RCP 6.0 in 2050s, RCP 4.5 and RCP 8.5 in 2070s whereas SH and MSH dominate HSH in five cases i.e. RCP 2.6, RCP 6.0 and RCP 8.5 in 2050, RCP 2.6, RCP 4.5, RCP 6.0. Some patches of SH will remain in Pak-Afghan border, Kashmir region and Himachal area in RCP 2.6 scenario in 2050s. In other scenarios, only MSH will exist in Pak-Afghan, Kashmir and Annapurna regions where some patches of HSH and HS will be available in Himachal areas except in RCP 2.6 in 2070s ( Supplementary Fig. S2).
When the entire Himalaya was examined collectively, we observed that the HSH and SH will increase in all scenarios (RCP 2.6, RCP 6.0, RCP 4.5 and RCP 8.5) during 2050s whereas HSH and SH will increase than now in two scenarios (RCP 4.5 and RCP 8.5) during 2070s and decrease in remaining scenarios (RCP 2.6 and RCP 8.5) (Fig. 5). Except in 8.5 scenario in 2070s, HSH and SH will diminish than 2050s. MSH will maintain similar pattern between now and 2070s in all scenarios whereas MSH will increase in three scenarios (RCP 2.6, RCP 6.0 and RCP 8.5) in 2050s. Area of MSH will decrease in 2070s than 2050s in all scenarios except RCP 4.5. More habitat is expected to be available for KMD in Uttarakhand, west Nepal and their associated areas in Tibet in future. Both HSH and SH will dominate Uttarakhand, west Nepal and Tibet in 2050s and 2070s. protected areas and refugia. A small fraction of KMDs range now or likely in the future (Fig. 6) are being protected by the currently designated network of protected areas (green shade in Fig. 7). Currently, only 17% of range inside the protected areas contains habitat with 0.5 or greater suitability. Many of the protected areas are located north of KMD's currently suitable habitat. As suitable habitats shift northward, more of the future range will fall outside currently protected areas (12-16%). If a site provides current and predicted future suitability, it carries higher conservation value because it could serve the purpose of a climate refugia potentially allowing that local population of KMD from needing to migrate to survive. By overlapping predicted future species range maps with the current ranges, we show that 22% and 24% of current species range will serve as climate refugia in 2050 and 2070, respectively (Fig. 7). Unfortunately, only 4% of these potential refugia are located inside present protected areas.

Discussion
Distribution of KMD: A game of moistures. The growth limiting factor for the species of trees in KMD (Kashmir Musk Deer) habitat is the moisture available during the pre-monsoon season [43][44][45][46][47] . In the southern Himalaya, moisture comes from precipitation arising from the Indian Ocean mainly during summer, making high rainfall in the eastern Himalaya leaving the western Himalaya drier than eastern Himalaya. During winter, the westerly winds from the Mediterranean sea brings more precipitation in the form of snow to the western Himalaya, but less to the eastern Himalaya 48 . The massive mountain range acts as a barrier to air flow from south to north, resulting in a dry and treeless Tibet 49 . Most of the coniferous trees respond negatively with temperature when temperature increase above 6.6 °C and positively with precipitation during pre-monsoon [50][51][52] . Spruce Picea smithiana, fir Abies pindrow, yew Taxus wallichiana, bluepine Pinus wallichiana, Himalayan birch Betula utilis and rhododendron Rhododendron campanulatum etc. are commonly recorded tree in the habitat of KMD. Growth of these trees are limited by moisture in pre-monsoon [43][44][45]53 . Our finding also suggests that Himalayan trees are more responsive to changes in the regime of precipitation than temperature. The strong influence of precipitation to KMD distribution likely results from the sensitivity of forests supporting KMD sites to adequate rainfall. The distributional boundaries we predicted are echoed by recent studies: Singh et al. 35 indicated Mustang, central Nepal as the eastern boundary and Ostrowski et al. 23 indicated Nuristan, Afghanistan as the western boundary. However, the distributional range of KMD as depicted by IUCN is markedly different from our data-driven modeling studies IUCN 32,35 . future distribution of KMD. West Nepal, Uttarakhand, Kashmir region, and Pak-Afghan border will support highly suitable and suitable habitat. Most of the current habitat of KMD will disappear in the 2050s or the 2070s in all climatic scenarios (RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5) except Uttarakhand and west Nepal www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ and their adjacent areas in Tibet. However, few patches of suitable habitat will remain in Himachal and close to Pak-Afghan border in 2050s in RCP 2.6 scenario. Mean temperature is predicted to rise by 1-2 °C in 2050s and by 1-3 °C in 2070s in four RCP scenarios 54 . Although temperature will keep rising, precipitation is predicted to behave differently depending on the amount of moisture in the air, where it comes from, and the convergence of air currents 55,56 . In the range of KMD, sites that are unsuitable today to the north of Uttarakhand and western Nepal will receive more pre-monsoon moisture in future (see Supplementary Note). This explains why those currently uninhabited areas show a suitable sites in future climates (Fig. 4).
Himalayas in Uttarakhand and west Nepal are not higher as Himalayas in Annapurna region, Pakistan and Kashmir region. Mountains above 8000 m from mean sea level are located either in Pakistan or Annapurna region and east of it. Not surprisingly this huge area and its diverse topography manifests in a very complex climate system that is not well understood 57 . This distinct dynamism of temperature and precipitation due to geographical differences between the western regions (Kashmir and Pak-Afghan) and the eastern region (Uttarakhand and west Nepal) are enough to justify the reason behind loss of habitat of KMD west of Uttarakhand (Kashmir region, Afghanistan and Pakistan) in future (2050s and 2070s).
Contrary to our finding, Khadka et al., 2017 found temperature play a greater role in distribution of musk deer. Therefore, Khadka et al. 58 predicted (i) that vast stretches of currently suitable habitat of Himalayan musk deer exist at very high elevation including above tree line and snowline although we know that the species needs a forest to survive, and (2) that this suitable habitat expands dramatically in future to higher elevation areas without any contraction. The unrealistic future projection of musk deer by Khadka et al., 2017 has another layer of complexity: plants do not respond to newly available suitable climate above treeline and snowline because soil may not have formed in enough quantity to support the plant growth in such large area 59,60 . Collectively, these lines of arguments and evidences support our finding that current species range of KMD will respond to future climate in a complex fashion, which is different from a species that is primarily controlled by temperature. We show that all kinds of habitats in Kashmir region, Himachal and Pak-Afghan border will be the ones showing the first sign of climate change; they will decline in the 2050s and 2070s.
There are a few reasons why prior modeling studies (e.g., Khadka   www.nature.com/scientificreports www.nature.com/scientificreports/ species and it changes with season. Different species can look alike as well as individuals of the same might look different 35 . Fecal pellets are not uniquely identified at a species level. Occurrences based on opportunistic sightings and anecdotal information has the potential to hamper the model as such occurrences have a higher chance of being incorrect or falling in lower quality habitat than our systematic survey of latrines. Khadka et al., 2017 and Khadka and James 2017 concluded that temperature (relative importance = 74%) is the most important variable influencing Himalayan musk deer distribution. In water scarce environments, a small difference in moisture availability can make a difference in habitat quality. Such ecological expectations can also be inferred from the fact that forest cover is an important factor to KMD given their timid nature 30 and vegetation growth is largely influenced by moisture in dry environments [60][61][62][63][64] . These deer occur sparsely if there is inadequate dense vegetation cover (>40% crown cover) to hide within during the day 30 . protected, non-protected habitats and climate refugia. We observed that only 17% of the current suitable range (projected probability >0.5) falls inside protected areas (PAs). In the future, this will decrease to 12-16%. The range are in low fraction and other researchers have reported similarly low fraction of range protected inside conservation area for other species of Himalayan range such as 19% for snow leopard 61 . Large portions of the suitable habitat of musk deer are located outside PAs where anthropogenic pressure are highly www.nature.com/scientificreports www.nature.com/scientificreports/ threatening biodiversity 65,66 . Only 4% of habitat in 2050s and 2% of habitat in 2070s will act as climate refugia inside PAs. However, 25% of current suitable range in 2050s and 16% in 2070s will act as climate refugia in entire range of KMD. climate refugia conservation; next generation conservation. Musk deer evolved 4 million years ago in Tibetan plateau 67 and survived until now despite various climate changes during the Quaternary period as well as after remarkable geological process in Asia including upliftment of Himalaya 68,69 . Climatic refugia are a plausible explanation for its survival through different climate changes, allowing the species later to disperse in the mountains of Asia when situation slowly normalized. Another plausible explanation for the survival relates to the ability of a species to cope with the changing environment and its interaction with other species. Unfortunately, the velocity of the climate change is faster than the response of the species is 70 . Therefore, the role of climate refugia is always crucial. The extinction of megafauna, such as woolly mammoth Mammuthus primigenius and Giant deer Megaloceros giganteus during the late Quaternary period was mainly because they could not find such refugia to survive through the glacial-interglacial cycle 71 . In the case of KMD, about quarter of the current range can serve as climate refugia in future (2050s and 2070s). These refugia are situated in Uttarakhand, west Nepal and west of Annapurna region.

Material and Methods
Study area: the greater himalaya. The greater Himalaya are often called the roof of the world and they extend from northern Pakistan to north east India throughout Nepal, Bhutan and the southern part of Tibet. Western disturbances from the west arise from the Mediterranean Sea and dominate the climate of western Himalaya during winter. During summer, the dominant weather pattern is driven by moisture arising from the Indian Ocean 48,72,73 . In Uttarakhand and west Nepal, the climate is different than Kashmir and Pak-Afghan as it is influenced by both south west monsoon and westerly storms 74 . The different weather phenomena in western and eastern Himalaya has resulted in a drier west and wetter east 75  presence records of musk deer. Following the genetic study of the samples collected in the Nuristan region of Afghanistan and the Mustang region of Annapurna Himalaya range in Nepal (Singh et al., 2019), Nuristan is confirmed as the western limit and Mustang as as the eastern limit of the KMD range. We determined the geographic coordinates of latrines between these two distribution limits of KMD in South Asia. We documented 136 occurrence records of KMD, 123 of them being in our primary collection area within five protected areas in Nepal (Annapurna Conservation Area, Dhorpatan Hunting Reserve, Shey Phoksundo National Park; Rara National Park, Api Nampa Conservation area), Musk Deer National Park in Pakistan, Nuristan of Afghanistan, and four protected areas in India (Dachigam National Park, Kedarnath Wildlife Sanctuary, Nanda Devi Biosphere, Govinda Pashu Vihar Wildlife Sanctuary). In total, we collected 45 latrine samples in Nepal, 61 in India, 21 in Pakistan, and 9 in Afghanistan (Fig. 1).
Latrines, which serves as an indisputable evidence of deer presence, can easily be distinguished from the excreta of other animals; they have the heaps of old and fresh pellets with a musky smell 42,[79][80][81] , and the pellets are much smaller and cylindrical compared to that of goat Capra aegagrus and sheep Ovis aries. Musk deer are shy, nocturnal and crepuscular forest-dweller species 82 . During the day, they hide in vegetation. Hence, musk deer are mostly inconspicuous to human. However, musk deer develop latrines by defecating repeatedly at a single site to maintain communication for various proposes 81,83,84 .
Geographical locations of latrine sites were determined first in Mustang, Nepal. In these expeditions, we performed systematic surveys for understanding habitat by laying out quadrats of size 20 m × 20 m along transects at different elevations. In the summer and autumn of 2016-2017, we visited several areas up to 4500 m in the high Himalayans of Nepal (Manang, Mustang, Kaski, Bhimthang of Annapurna region). Occurrence records from remaining four protected areas of Nepal were collected by random sightings of latrine sites and animals by the game scouts and rangers working in each protected area. For India and Pakistan, we used a retrospective study design. Coordinates of latrine sites systematically collected for Alpine musk deer from Uttarakhand, India and from Pakistan by co-author of this paper from 2003 to 2017 were considered to be KMD instead of Alpine musk deer. Data in Pakistan were collected from systematically distributed quadrats. Additionally, published literature was reviewed and presence locations of musk deer from Afghanistan and Dachigam National Park, Kashmir region, India were extracted by scanning and georeferencing previously produced range maps (Fig. 1). The data obtained from previously published sources (secondary data) comprised of less than 10% of our pool of occurrence records whereas our primary data collection effort contributed over 90% of the occurrence records. potential predictors. The current climate is represented by raster layers for 19 climatic variables with the resolution of 30 arc seconds (approximately 1 km) ( Table 1)  www.nature.com/scientificreports www.nature.com/scientificreports/ KMD, we supplemented this standard set of potential covariates with additional relevant variables (Table 1). In total, 26 covariates were applied to the first set of models.
Our initial modeling exercise indicated that several of these covariates secured a very low variable importance. Some of these minimally important covariates also had other problems. Specifically, "distance to the nearest source of water" required computing distance from each of our occurrences and background points to the water bodies and identifying the nearest source. This turned out to be computationally intensive. So, we safely dropped minimally important covariates "distance to the nearest source of water". Further, we dropped "aspect" "altitude" and "land cover" because of their static nature and they could lead to over fitting the model. We also adapted filter i.e. approach to eliminate over fitting. This approach was useful in optimizing the model coefficients such that a more stable model is trained. Specifically, we tested spatial transferability of the models by comparing various model evaluation metrics [86][87][88] . The goal here is to find a model that performs in both training and testing regions similarly. With various regularization multipliers, we found that a value of 2 minimizes the delta AIC (difference in AIC of the models in training vs testing regions). Finally, we selected 22 variables for processing our model to understand present distribution of the KMD. To forecast range patterns in the years 2050 or 2070, we used the same variables examined the present distribution of KMD in this study and applied the four climate change scenarios as used in Fifth Assessment by the International Panel on Climate Change (Representative Concentration Pathways RCP2.6, RCP4.5, RCP6, and RCP8.5). For forecasting future distribution, we used Global Circulation Model (GCM). We reviewed various literatures related to GCM. BCC-CSM1-1, CCSM4, HadGEM2-CC and MIROC 5 have been mostly used in SDM to perform future prediction. All of these models are also used in Fifth Assessment IPCC report (IPCC AR5). CMIP5 publications database reveals that CCSM4 is at top level which has been used in 440 articles, which is followed by MICROC 5 (409) and then by BCC-CCM1-1 (389). HadGEM2-CC has been used by 337 publications. We selected BCC-CSM1-1 because it is based on CCSM and itself widely used for projecting future distribution 89 . Therefore, using BCC-CSM1-1 provides double benefit over other GCMs. Other reasons behind applying BCC-CSM are i.) BCC-CSM1-1 is developed to issue global climate predictions and impact assessments at monthly, seasonal and inter-annual time scales, particularly over East Asiaand ii.) to do research on climate and climate change issue 90 . Our future predictions are related to climate change issues. We found BCC-CSM1-1 most fit in our case. All variables were resampled with the resolution of 30 arc seconds to use in model.

Maximum entropy (Maxent). MaxEnt is a widely used Species Distribution Model (SDM) for predicting
species distribution in the fields of ecology and conservation 64,[91][92][93] . MaxEnt 91 operates by computing the ratio of the probability density function of multivariate environmental space of covariates in observed locations to that of the entire study area and it is particularly suitable for modeling presence-only data. It is mathematically related to Poisson Point Process 94 , another model suitable for presence-only analysis. Elith and Leathwick (2009) consider MaxEnt to be a particularly efficient model and we determined it was a good choice for modeling the distribution of KMD.
An important assumption of SDM is that distributional data of a species can be collected by either randomly or systematic sampling of occurrence records 95 . Our occurrence data were collected systematically in Annapurna Conservation Area (Annapurna region), Api Namba Conservation Area in Nepal, Nanda Devi Biosphere, Govinda Pashu Vihar Wildlife Sanctuary but the same scheme was not used in Afghanistan, Kashmir region and far western Nepal. Therefore, it can be argued that the occurrences might have some level of sampling bias. When background points are drawn with similar bias as presences, this, in theory, dampens the impact of sampling bias in presences. As expected, Phillips et al. 2009 showed improved model performance with biased background drawing. The bias can be directly obtained from the intensity of presences in the grid cells but making sure that the grid cells without any occurrences still contribute the background points with a much smaller probability 96 .
We sampled background points from occupied grid cells with a bias equal to the intensity of presences. For unoccupied grid cells, we sampled background points at two probabilities: 0.01 and 0.1, and compared the predicted surface and evaluation metrics of these three schemes of background drawing (two discussed above plus one without applying any bias in the entire study area). We checked current range of the species resulted from three types of biases with the precisely known habitat of the species in reference to Himalaya range. We found that the results slightly varied among the three types of biases (AUC = 0.99 for 0.1 bias, 0.96 for 0.01 bias and 0.99 without bias). We will report results only for 0.1 bias because this secured highest AUC score. Following Guo et al., 2017, we classified the continuous probability into four categories: unsuitable (0-0.2), marginally suitable (0.2-0.5)(MH), suitable (0.5-0.7)(SH), and highly suitable (0.7-1.0)(HSH) 60 to understand current and future distribution of KMD.

Data availability
Musk deer is endangered species and poaching for the musk is major cause of population decline of the deer. Our presence records are location of latrine sites and it was found that poachers rely on freshness of latrines site to track the deer. These sites are used repeatedly by musk deer for many years. Therefore, we can not make our presence data publicly available. However, we can provide these data upon the formal request from the researchers.