Using an ensemble modelling approach to predict the potential distribution of Himalayan gray goral (Naemorhedus goral bedfordi) in Pakistan

Journal Pre-proof Using an ensemble modelling approach to predict the potential distribution of Himalayan gray goral (Naemorhedus goral bedfordi) in Pakistan Shahid Ahmad, Li Yang, Tauheed Ullah Khan, Kunyuan Wanghe, Miaomiao Li, Xiaofeng Luan PII:

S2351-9894(19)30350-6

DOI:

https://doi.org/10.1016/j.gecco.2019.e00845

Reference:

GECCO 845

To appear in:

Global Ecology and Conservation

Received Date: 22 June 2019 Revised Date:

29 October 2019

Accepted Date: 9 November 2019

Please cite this article as: Ahmad, S., Yang, L., Khan, T.U., Wanghe, K., Li, M., Luan, X., Using an ensemble modelling approach to predict the potential distribution of Himalayan gray goral (Naemorhedus goral bedfordi) in Pakistan, Global Ecology and Conservation (2019), doi: https:// doi.org/10.1016/j.gecco.2019.e00845. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Using an ensemble modelling approach to predict the potential distribution of Himalayan gray goral (Naemorhedus goral bedfordi) in Pakistan Shahid Ahmad1✝, Li Yang1✝, Tauheed Ullah Khan1, Kunyuan Wanghe1 Miaomiao Li1, Xiaofeng Luan1a,* 1 School of Nature Conservation, Beijing Forestry University, NO.35 Tsinghua East Road Haidian District, Beijing, 100083, P. R. China ✝ Joint first authors/these authors contributed equally to this work * Correspondent: Xiaofeng Luan School of Nature Conservation, Beijing Forestry University, NO.35 Tsinghua East Road Haidian District, Beijing, 100083, P. R. China. (E-mail: [email protected]; Tel: +86-13910090393; Fax: +0086-10-62336724) Keyword: Niche Change; climate change;;;; Habitat change; Habitat prediction; ungulates; Capture-mark recapture

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

Abstract Global warming has negative impacts on the distribution of large ungulates, particularly for species occupying narrow distributional ranges. Knowledge of how climate change will affect future distributions is imperative for designing effective conservation action plans for at risk species such as the Himalayan gray goral (HGG), a cliff-dwelling mountainous goat. We sought to evaluate the potential distribution of Himalayan gray goral (HGG) under future climate change scenarios using ensemble modeling approaches. HGG data were obtained from previous published surveys, publications, and occurrence records ranging from 1985 to 2018. we also conducted survey in 2017-2018 using double observer method based on capture mark recapture. (Suryawanshi et al., 2012; Tumursukh et al., 2015). Later on we double check the record and remove double observation. After quality control screening, 139 records remained for analysis. Resulting species distribution models (SDMs) results showed sufficient internal evaluation metrics, with all TSS values being > 0.7. The random forest (RF) modelling technique had on average the lowest true skill statistics (TSS) value, However the multivariate adaptive regression splines (MARS) modelling technique had the highest. The ensemble modelling internal evaluation metrics indicated adequate results with values ranging from 0.827 to 0.843. Annual mean temperature (Bio1) and annual precipitation (Bio12) were found to be the most important climatic variables impacting the potential distribution of HGG. HGG habitat determined to be suitable in both current and future climate scenarios decreased in all Representative Concentration Pathways (RCPs) scenarios with the exception of RCPs 2.6. Suitable habitat in both current and future climate scenarios remained consistent in the time periods of 2050 and 2070 under RCP4.5 while fluctuating in 2030, 2050 and 2070 under RCP 2.6. However,, the suitable habitat under current and future scenarios declined in 2030 under RCPs 4.5 and in 2030, 2050, and 2070 in scenarios RCPs 8.5. Currently suitable HGG habitat was located in an area where the species is known to be locally extinct. Further work is necessary to determine the key drivers of local extinction events in an effort to mitigate population crashes. Our work will assist in formulating conservation actions for the HGG in the context of climate change, and provide a platform for continued monitoring efforts of the species.

29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68

Introduction The ecological niche of a species is the interaction of space and condition that promotes viable populations and species perseverance through time (Choudhury et al., 2016). Factors associated with occupied habitat niches explain species origination, persistence, distribution and competition with other sympatric wildlife for space and resources (Silvertown, 2004). Temperature, precipitation and biological interactions are principal drivers known to impact species distribution (Woodward, 1987; Abolmaali et al., 2018). Such drivers can be largely impacted by global warming, which may cause population crashes that result in an endangered status or extinction (Parmesan, 2006; Loarie et al., 2008; Cuena-Lombraña et al., 2018). Climate change effects have been documented for almost all ecosystem types and species, including human (Chang et al. 2015; Iannella et al. 2018; Ruete & Leynaud 2015) animals (Iannella et al. 2018; Ruete & Leynaud 2015), plants (Ashraf et al. 2018; Wei et al. 2018). The situation is expected to be worsed as the speed and magnitude of global environmental warming is unceasingly accelerating (Meehl et al. 2007). Developing countries are particularly at risk for negative outcomes associated with climate change despite their more negligible contributions to global warming because. Pakistan is one such developing country facing the threatening consequences as a result of climate change, particularly the more remote and mountainous regions of the country occupied by the high-elevation Karakoram, Hindu-Kush, and Himalayan mountain ranges. Himalayan ecosystems are some of the most sensitive to global warming in the world (Laghari JR 2013). They harbor a number of unique flora and fauna, including a number of mountainous ungulates. Unfortunately, many such species lack adaptive strategies to novel climatic conditions (Huntley et al. 1995; Pearson & Dawson 2003; Yates et al. 2010). Knowledge of potentially suitable areas where large ungulates occupying a narrow niche may find refuge under future climatic conditions can enhance conservation effectiveness and reintroduction efforts (Dubuis et al., 2011; Kaky and Gilbert, 2016). Such conservation action plans are appropriate for the Himalayan gray goral (Naemorhedus goral bedfordi), a cliffdwelling rupicaprine, classified as Near Threatened throughout its range by the IUCN Red List (Abbas et al., 2015). HGG inhabit rugged high-altitude terrain at elevations ranging from 649 to 3789 m. The species distributed from Pakistan southwestward to Nepal. (Ashraf et al. 2017; Roberts & Bernhard 1977). HGG serve as a criticaly important prey item for a number of other carnivore species, including snow leopard (Panthera uncia), common leopard (P. pardus), and Eurasian lynx (Lynx lynx). Despite its importance, the HGG is grossly understudied, with very limited understanding surrounding habitat selection and potential distributional shifts. Data analyses and statistical tools, such as species distribution models (SDMs) can be used to address such knowledge gaps. Of the studies outside of Pakistan, 37.1% focused on the distribution and estimation of population abundance. When considering the dataset in its entirety, only 2.7% of articles addressed habitat suitability and climate change impact on HGG. (See detail in Appendix I).

69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106

SDMs are important tools that combine presence and pseudo-absence data with abiotic data (Elith & Leathwick 2009a; Thuiller et al. 2009). SDMs are used to predict the geographic distribution of species for ecological, biological and resource management (Ashraf et al. 2018). Discrepancies between differing SDMs make the selection of the most appropriate model difficult (Elith & Leathwick 2009a; Elith et al. 2011; Renner & Warton 2013). This is particularly true when models are used to project the distribution of species into independent situations, such as the projection of species distribution into future climatic scenarios or geographic space (Thuiller 2004; Thuiller et al. 2009). The ensemble modeling (ENM) approach is one such method to circumvent these issues. ENM comprehensively evaluates species distribution and its potential changes under climate change (Thuiller et al. 2009), with BIOMOD2 serving as a suitable platform (Cianfrani et al. 2011; Thuiller et al. 2009). To predict the impacts of climate change on species distribution, multiple scenarios must be considered (Araújo & Rahbek 2006; Beaumont et al. 2008; Beaumont et al. 2007; Bellard et al. 2012; Parmesan 2006; Sala et al. 2000) based on various Representative Concentration Pathways (RCPs, IPCC, 2013). RCPs explain scenarios based on an assumption of socio-economic condition, greenhouse gases, and air pollutant emission to provide trajectories for climate change (Albuquerque et al. 2018). This study sought to apply an ENM approach to investigate (1) the impact of several climate change scenarios on future HGG habitat suitability, (2) the species ecological niche concerning climate change, and (3) the present and future distributional areas of HGG across Pakistan. While, numerous studies have examined the population distribution, (Abbas et al. 2015), habitat suitability (Roy et al. 1995), and feeding ecology of the species (Chaiyarat et al. 1999), this is the first study to our knowledge concerning the geographical niche, habitat suitability, and changes in spatial geographic distribution under climate change scenarios for HGG.

107

Data collection

Materials & Methods Study area The study area (33-360 N and 71-730 E) encompassed the entirety of the known HGG distribution in Pakistan with elevations ranging from 649 m to 3789 m (Appendix II). Climatic conditions are highly diversified, oscillating between monsoon- influence moist temperate zone and semi-arid cold-deserts (Jacobose 1993). Four vegetation zones correspond to altitudinal gradient - snowfields, alpine meadows, subalpine scrub zones and dry alpine steppes (Bischof et al. 2014). The area house a number of large mammals, including Asiatic ibex (Capra ibex sibirica), markhor (Capra falconeri), snow leopard ( (Panthera uncia) , lynx (Lynx lynx), common leopard (Panthera pardus), grey wolf (Canis lupus), wild boar (Sus scrofa) and brown bear (Ursus arctos) (Bischof et al. 2014; Virk et al. 2003). Figure 1.

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147

Occurrence data of HGG were compiled from previous literature published from 1985 to 2018 (Drake & Beier 2014; Wei et al. 2018; Xu et al. 2008; Yang et al. 2016; Zhang et al. 2016). For the analysis, only data with precise location information were selected (Wei et al. 2018). Records that conflicted with the dataset, were information deficient, or lacked detailed descriptions were excluded. Information from duplicate locations were removed to reduce geographical bias and spatial autocorrelation (Beck et al. 2014; Marks 2011; Ponder et al. 2001). The ‘thin()’ function using the “spThin” package in R software version 3.5 was used with 100 iterations for a given thinning distance of 1 km within HHG home range in Pakistan to reduce bias caused by occurring agglomerate points. Locality information was then extracted from the documented records. Valid records were assigned with coordinates from public resources such as Google Earth v 7.1.2. All valid records were cross-checked by habitat information using Google Earth v 7.1.2 and Land-cover (GLCNMO, version3) (Tateishi et al. 2014). Finally, 139 records were selected for modeling. Pseudo-absence data can be important for ENM, such dataset associate to collection bias (Ponder et al. 2001). Therefore all the occurrence points in the study area from 1950 onward were downloaded for all species, except HGG, using the GBIF Global biodiversity information facility (GBIF https://www.gbif.org) in the “rgbif” package in R software version 3.5). A total of 4,000 points were randomly selected to generate pseudo-absence data. Capture Mark Recapture: Capture-mark recapture (CMR) surveys using the double observer’s survey method (Suryawanshi et al., 2012; Tumursukh et al., 2015) were completed in areas of Pakistan where HGG were known to occupy but lacked any published literature on occurrence. Surveys were conducted from September 30 to October 15, 2017 and the first week of January 2018. Survey areas were divided into small blocks no more than the daily movement of the HGG in km2. High altitude ridgelines were considered boundaries where there was less chance of crossing of animals. Observers (A and B) both visually scanned blocks, with either Observer A 30 minutes in step ahead of Observer B, or Observer A adopting a different route than Observer B within the same block (Tumursukh et al., 2015). Study areas were scanned for HGG from 7.00 am to 10.00 am and then from 3.00 pm to 6.00 pm when the species is known to be foraging (Jackson et al., 1996; Schaller, 1977) using 8×25 binoculars Pentax (XCF) and spotting scopes 20×60 (Nikon). GPS (Garmin 62s) was used to record the coordinates. Upon sighting of a herd, individual HGG were divided based on sex and age class. Upon completion of the scanning session, Observers A and B compared their data and similar groups were identified on the basis of herd size, demographic categories, habitat types and location and similar and different groups sighted were identified. Then the similar data were discarded Environmental Variables. The Environmental variables used to characterize the ecological niche model were selected based on climate and topography. Multicollinearity may hinder relationship interpretation of predictor variables (Heikkinen et al., 2006). Principal component analysis (PCA) was completed in R (version 3.5) to identify and remove variables exhibiting significant multicollinearity. A total of

148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187

six climatic and topographic parameters remained - mean annual temperature (Bio1), mean diurnal range (Bio2), the maximum temperature of the warmest month (Bio5), annual precipitation (Bio12), precipitation of driest month (Bio14), annual heat moisture index (AHM), and distance to water. The three greenhouse gases (GHS) concentration trajectories ranged from low to high concentration (RCP2.6, RCP4.5 and RCP8.5 w/m2) to correspond to enhance in global radiative forcing values in the year 2100 relative to pre-industrial values. It was assumed that global mean temperature would increase by 0.30°C to 4.8 0°C by 2100 across all RCPs (Van Vuuren et al. 2011) Climate data, including precipitation, mean temperature, minimum temperature, and maximum temperature, were obtained from CCAFS (http://www.ccafsclimate.org) and ESGF (Coupled Model Intercomparison Project phase5 https://pcmdi9.llnl.gov/projects/esgf-llnl). Environmental variables for 12 months per 20 years 2030 (2020 - 2040), 2050 (2040 - 2060) and 2070 (2060-2080) from 1985 to 2018 (obtained from CRU TS v. 4.01) were calculated by ANUSPLIN v4.36. This package provides the thinplate smoothing spline model for creating climate layers. The R “dismo” package was used to derive the nineteen bioclimatic variables for month wise temperature and precipitation. Topographical data was taken from SRTM digital elevation data (http://srtm.csi.cgiar.org/). Distance to water was created by Global Surface Water (Hansen et al. 2013) following Wang et al., (2012). The annual heat and moisture index was calculated using the formula by (BIO1 + 10)/ (BIO12/1000). All variables were obtained with the same extent and spatial resolution of 30 sec (cell size approximately 1 km2) (Appendix III). (Araújo & Rahbek 2006; Dussault et al. 2006; Franzmann 1981; Piao 1995; Wattles & DeStefano 2013) Table 1. NOTE: All variables were tested by variance inflation factor (VIF) (threshold = 10; any variable has VIF value higher than 10 was excluded).. AHM were calculated by (BIO1 + 10)/(BIO12/1000). Species Distribution Modeling. SDMs are commonly used in the modeling of species distribution under climate change. (Elith et al. 2006; Elith & Leathwick 2009b; Guisan & Thuiller 2005; Guisan & Zimmermann 2000; Hijmans & Elith 2013; Miller 2010). Recent advancements, such as the use of computer based algorithms, in the calibration of SDMs (Elith et al. 2006; Wisz et al. 2008) in conjunction with varying ensemble procedures (Araújo & New 2007; Buisson et al. 2010; Kindt 2018; Luedeling et al. 2014; Marmion et al. 2009; Thuiller et al. 2009) allow for more robust analyses that produce more reliable habitat suitability maps. Ten modeling techniques were used: Random forest (RF), generalized boosting modeling (GBM), generalized linear modeling (GLM), artificial neural network (ANN), classification tree analysis (CTA), surface range envelope (SRE), flexible discriminate analysis (FDA), maximum entropy (MAXENT.Phillips), multivariate adaptive regression splines (MARS), and generalized additive models (GAM). The biomod model parameters were modified to determine the best performance technology for further research. Each SDM (with each modeling algorithm, split the data into 80% and 20% for training and testing, respectively) was evaluated using the true skill statistic (TSS) (Allouche et

188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219

al. 2006) with ten evaluation runs. ROC and TSS were used as two modes of criteria for model effectiveness. RF, MARS, GBM, and MAXENT.Phillip were selected for ENM (Fig 1). Model results were selected by threshold (mean TSS value among all evaluation runs) for ensemble modeling to reduce uncertainty. Potential distribution was calculated from an ensemble set of models or predictions, and a result of mean probabilities across predictions was selected. Then, the potential distribution under climate change for the three periods was projected. We repeated this process three times, resulting in 30 SDMs in total (3 ensemble results for each period in each RCP) under current and future climate conditions. For each period under each RCP, the final SDM result was calculated by the weight as determined by TSS. The cut-off value calculated by BIOMOD2 was considered as a threshold to transform the models from environmental suitability into presence: output value higher than the threshold was considered “present”; otherwise, “absent.” Habitat Change. We overlayed potential habitat under the current climate condition with each period under each RCP. Areas which served as suitable habitat in both current and future climate scenarios (each period under each RCP) were considered as “always suitable areas”. Areas that were suitable in current climate scenarios but lost in future climate scenarios were categorized as “habitat loss area” and areas that were unsuitable in current climate scenarios but gained in future climate scenarios were categorized as “habitat gain areas”. Spatial analyses were carried out in ArcGIS (ver. 10.2; ESRI Inc., Redlands, CA, USA) and Excel 2013 (Microsoft Corp., Redmond, WA, USA) Results. Model Performance. The SDM results showed proper internal evaluation on average. All TSS values were greater than 0.72 (Appendix III Table S2.). MARS had the highest average TSS values amongst the other Algorithm. To increase model accuracy, the average of all TSS values was set as the threshold for building ensemble models. Model selection was based on ROC and TSS values (ROC > 0.92 and TSS ≥ 0.72) (Fig 2). The response curves of a model independently of the algorithm used for GBM, MARS, MAXENT. Phillips and RF are presented in Appendix III (figure S3~figure S6). This showed that Bio1 (annual precipitation) and bio12 (mean annual temperature) are the main factors which highly affecting the HGG distribution. Figure 2.

220

Temporal Niche Changes.

221 222 223 224 225 226

Two main components resulted from the PCA. PC1 explained 50.4% variance of the current climate scenario while PC2 explain 30.16%. The two axes of PCA explain the distribution of climatic variable distribution and percentage of inertia. Black dots represent the current situation in each RCPs. Red dots represent HGG distribution under bio1, bio12 in 2030 in each RCP. Green dots showed HGG distribution under the influence of bio1, bio12 in 2050 in each RCP and Blue dots revealed HGG distribution under bio1, bio12 climatic variable in 2070 in each

227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266

RCP. With two main gradients, PC1 is mainly corresponding to the temperature and PC2 to precipitation. The correlation circle revealed the variance contribution of each PC1 and PC2 for each RCP. The overlap between the partial and the total niche also tends to be stable for different RCPs and different periods (Fig 2). The results showed that Niche was statistically significant for all comparison (current to future). Figure 3. Variable Importance Measures. Bio1 (annual mean temperature) and Bio12 (Annual Precipitation) were the primary variables found to influence HGG habitat suitability, 64±0.05%, and 36±.09%. Each variable relative influence is given in Table S2 (Appendix III). Among the six selected environmental variables, annual mean temperature (Bio1) and annual precipitation (Bio12) had the greatest effect on HGG distribution under the examined climate change scenarios followed by max temperature of warmest month (Bio5), precipitation of driest month (Bio14), Annual heat moisture index (AHM) and water distance in different RCPs. Precipitation of driest month (Bio14) had the highest variability in the future climate scenario. For complete detail of environmental variable variability in each RCP in each time period see Figure S2 variable change in the near future (Appendix III).Predicted Future Suitable HGG Habitat under Global Warming Scenarios The Earth’s changing climate affects species range and suitable habitat. Ranges and suitable habitat can move, shrink, or grow as a result of climate changes. Suitable HGG habitat fluctuated based on the RCPs and time frame examined. However, the overall range of suitable HGG habitat shrank in future climate scenarios compared to the current suitable habitat and range of HGG. The predictions for 2030, 2050 and 2070 (under the Rcp2.6 Rcp45 and Rcp85) climatic models showed almost or slightly less than the current scenario, with slight fluctuations depending on RCPs and time period examined. For the Rcp8.5 scenario, 21,533km2 of always suitable habitat was predicted in 2070, exhibiting the greatest loss of either parameter amongst all RCPs and time periods. Currently, the species is locally extinct in the southern mountainous area especially (District Mardan, Buner, Malakand and Nowshehra) while the area is suitable. It may be due to human activities and illegal hunting. It may also be results of climate changes impacts on flora distribution, which HGG used in different life span periods. Our study revealed that combining data, get from multiple sources can offer an opportunity to identify suitable habitat for targeted species on a local scale. Fig 2 showed the HGG complete detail of habitat and range gain and loss in the future in a different periods and in different climatic scenarios. Figure 4 Potential Habitat Change Current species range in the study area is 34,133 km2, with the elevation range of 649m to 3789 m. HGG range in all RCPs declined by 12.10%, 17.66%, and 25.33%, respectively. After projection into future and compared to the current distribution, the total area of the always suitable habitat for HGG under the RCPs (2030RCP2.6, 2050RCP2.6, 2070RCP2.6, 2030RCP4.5, 2050RCP4.5, 2070RCP4.5 and 2030RCP8.5, 2050RCP8.5, 2070RCP8.5) would

267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298

decline by 25%, 23.4%, 18% and 26%, 31.3%, 31.3%, 26%, 40%, and 57%, respectively. The always suitable habitat decline in all RCPs scenarios except in RCPs 2.6. Fig4 showed the potential distribution habitat change of HGG under climate changes in the future in different climatic scenarios. It showed that species is under pressure of both climate and human. As already discussed that in some areas, species is locally extent but the habitat is still suitable and vise versa. Figure 5. Discussion. Currently suitable HGG habitat primarily aggregated in the lower northern region of Pakistan. These results align with our field observations and the known suitable habitat cited in the literature (Ashraf et al. 2018, Abbas et al., 2015) and unpublished data. Understanding species response to climate change is crucial for formulating effective conservation solutions (Quintero & Wiens 2013; Thomas et al. 2004; Warren et al. 2008). Numerous studies have already warned the direct effect of climate change on ecosystems across a wide range of taxa. (Muñoz-Mendoza et al. 2017; Quintero & Wiens 2013; Walther et al. 2002). The current study illustrates the first large-scale assessment of changes in distribution of HGG under climate change scenarios. Model results showed considerable predictive performance with sufficient TSS values, showing high predictive capacity (Domíguez-Vega et al. 2012; Wogan 2016). In Pakistan, HGG is one of the least studied mountainous species (Shackleton 1997). We covered the entire known distribution area of the species to predict both current and future distribution using ensemble modeling approaches using different modeling techniques and various greenhouse gas emission scenarios across several time periods. Currently suitable habitat remains same in periods of 2050 and 2070 under RCP4.5 and fluctuates in 2030, 2050 and 2070 under RCP 2.6. While both suitable habitat and range decline in 2030 under RCPs 4.5 and in scenarios RCPs 8.5 (2030, 2050, and 2070). Impact of climate change on the distribution of large mammals is not homogeneous across the species. Fluctuation in the predicted habitat under climate change varied among the species. Thuiller et al. (2006) results suggested significantly distinctive responses of African mammals to global climate change scenarios, where species lost almost all suitable habitats, while others were quite stable or even gained new habitat. Our results are aligned with literature. (Levinsky et al. 2007, BarbetMassin et al. 2009, Ogawa-Onishi et al. 2010) Luo et al., (2015) have predicted about 30–50% loss of distribution area for 22 ungulates in

299

Tibetan Plateau. In the current study, annual mean temperature and annual precipitation are

300

highly affecting the distribution of HGG. Kabir et al., (2017) also showed the effect of annual

301

precipitation on species distribution in the northern area of Pakistan.

302 303 304 305

Models’ results also displayed some topographic -climatically suitable sites within middle and lower northern parts of Pakistan where the species was not found during the field surveys. These recent local extinctions occurred in suitable sites where the species was previously present with very small population size;the reasons may be attributed to stochastic events or human

306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345

disturbances (overgrazing poaching and illegal hunting). Nonetheless, further surveys efforts are encouraged to investigate which factors had contributed to it. Our results showed that under current scenarios, the HGG distribution range was more influenced by bio1 (annual mean temperature) and bio12 (annual precipitation). Our results showed that for the years 2030, 2050 and 2070 under the three tested RCPs (2.6, 4.5 and 8.5) the geographic distribution of HGG would shrink under the future climatic conditions. The projected models showed habitat range shifts through the disappearance of this species. The reason for the range shift is that the climatic envelope (precipitation and temperature) of HGG will become less suitable for survival. In our study area, HGG is residing in the isolated pitches; therefore, it is more vulnerable to climate changes because the species having narrow ecological niches is more sensitive to climate change than the larger ecological niches. (Khanum et al., 2013;Abolmaali et al., 2018). The large mammal’s response to global warming climate mainly depends on the physiological characteristic and the feeding source in the suitable habitat (Zhao et al., 2017). HGG is selective in feeding; (Ashraf et al. 2018). The plant species are highly sensitive to climate change than animals. Although in our study area the research work to figure out the climate change impacts on plant species has not been conducted but in near countries plant species faced reduction in suitable habitat and range for many species such as Myristica dactyloides in India (Remya et al., 2015), Fritillaria cirrhosain in China (Zhao et al., 2017), Artemisia aucheri, A. sieberi and Daphne mucronata in Iran (Sanjerehei and Rundel, 2017; Abolmaali et al., 2018). Apart from climatic changes scenarios, anthropogenic activities can also trigger the distribution of species. Anthropogenic activities may result, over-exploitation of natural resources and over-grazing, causing habitat loss. In Pakistan, the high dependency of local people on the natural resource is causing habitat loss, and it is expected that actual habitat loss could be higher than the predicted loss. Due to the lack of reliable data, it is hard to estimates the relative impact of climatic change and human activates. The current research work is the first attempt to model HGG in its distributional range, describing the relationship between environmental factors and habitat suitability and predicting habitat suitability changes under climate change. Our results showed that habitat suitability of HGG is highly related to climate change which areconsistent with previous studies of climate change in the same ecosystems (Kabir et al. 2017). We add on to previous research on climate change in general, showing for the first time the potential changes in the habitat of HGG, under future climate changes scenarios. Our results provide a strong reason to carry out further work because the HGG distribution is poorly documented. For conservation plans and strategies, it is strongly recommended to declare the predicted suitable areas as a nature reserve that harbors more diversity with few risks both in current and predicted future climatic condition. Implications for Conservation Management. Effective control methods and rational management approaches based on broad considerations are necessary to determine the relationship between species and the environmental factors (as climate change, land use, and anthropogenic activities) that will most

346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383

contribute to their distribution (Pyke et al. 2008). For HGG, we propose that habitat regions most sensitive to climate change be identified as target priority areas for species perseverance (Pyke et al. 2008). Distribution modeling is and efficient and cost-effective tool that can provide an early warning system to identify the potential habitat under a series of climate scenarios (Warren et al. 2008). This provides the knowledge base neceesary for local governments and conservation agencies to select future sutatible areas as candidates for nature reserve establishment. For HGG, we suggest the middle and lower northern area of Pakistan and Kashmir as prime candidate areas for HGG protection. Secondly, many studies have focused on the species having socioeconomic importance for the local people trophy hunting like Asiatic ibex, markhor, and blue sheep, while, HGG is ignored in this regard by the local government as well international NGOs. The combination of local effort and governmental strategies can improve management and conservation plans. Model limitations: The current study has several limitations, it needs to be clarified (1) The relationship between density and occurrence frequency of species remains unclear (Carbone et al. 2001; Li et al. 2010; Liu et al. 2013) (2) Different survey methodologies come with varying depending on species, especially for large mammals (3) A microclimatic variation on different time periods may cause spatial changes in the potential distribution of species (here; we only create variable annuals or seasonal population changes (Broennimann et al. 2014; Sedgeley 2001). (4) Using high-resolution environmental variables may cause potential bias (Wang et al. 2012) (5) Modeling limitations (e.g., part of a species’ fundamental niche may not correspond to any combination of environmental variables; species interactions may cause bias) (Cruz-Cárdenas et al. 2014; Elith & Leathwick 2009b; Guo et al. 2015; Maiorano et al. 2013). However, occurrence records from different resources already contain potential information on habitat, which could reduce bias caused by variables representative of the model. In addition, conservation action mainly addresses on species richness, which can also reduce risk from modeling performance and potential bias for single species. Conclusion In the current study, we modeled the potential distribution of HGG in current and future climate scenarios across its distribution range in Pakistan. The principal environmental variables (Bio1 and Bio12)had the greatest influence on HGG distribution. Overall, suitable HGG habitat range will decline under future climate scenarios in 2050 and 20170. The current modeling framework has shown the suitable habitat in current climatic condition as well as in future under difference greenhouse gases scenarios Our study could be helpful to formulate and implement the conservation strategies for the HGG to mitigate the climate change effects and human disturbance on its distribution

Competing interests. We declare that we don’t have a competing interest Funding; Ministry of Science and Technology of the People’s Republic of China (Research and application of key techniques on endangered species conservation and prediction of forest fire and pests in response to climate change, 2013BAC09B00) Reference Abbas F-i, Mian A, Akhtar T, and Rooney TP. (2015). Status and Future Management of Grey Goral (Naemorhedus goral bedfordi) in Pakistan. Journal of Bioresource Management 2:2. Albuquerque F, Benito B, Rodriguez MÁM, and Gray C. (2018) Potential changes in the distribution of Carnegiea gigantea under future scenarios. PeerJ 6:e5623. Allouche O, Tsoar A, and Kadmon R. (2006) Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of applied ecology 43:1223-1232. Araújo MB, and New M. (2007) Ensemble forecasting of species distributions. Trends in ecology & evolution 22:42-47. Araújo MB, and Rahbek C. (2006) How does climate change affect biodiversity? Science 313:1396-1397. Ashraf N, Anwar M, Oli MK, Pine WE, Sarwar M, Hussain I, and Awan MS. (2017) Seasonal variation in the diet of the grey goral (Naemorhedus goral) in Machiara National Park (MNP), Azad Jammu and Kashmir, Pakistan. Mammalia 81:235-244. Ashraf U, Chaudhry MN, Ahmad SR, Ashraf I, Arslan M, Noor H, and Jabbar M. (2018) Impacts of climate change on Capparis spinosa L. based on ecological niche modeling. PeerJ 6:e5792. Beaumont LJ, Hughes L, and Pitman A. (2008) Why is the choice of future climate scenarios for species distribution modelling important? Ecology letters 11:1135-1146. Beaumont LJ, Pitman A, Poulsen M, and Hughes L. (2007) Where will species go? Incorporating new advances in climate modelling into projections of species distributions. Global change biology 13:1368-1385. Beck J, Böller M, Erhardt A, and Schwanghart W. (2014) Spatial bias in the GBIF database and its effect on modeling species' geographic distributions. Ecological Informatics 19:10-15. Bellard C, Bertelsmeier C, Leadley P, Thuiller W, and Courchamp F. (2012) Impacts of climate change on the future of biodiversity. Ecology letters 15:365-377. Bischof R, Ali H, Kabir M, Hameed S, and Nawaz MA. (2014) Being the underdog: an elusive small carnivore uses space with prey and time without enemies. Journal of Zoology 293:40-48. Broennimann O, Petitpierre B, Randin C, Engler R, Breiner F, D'Amen M, Pellissier L, Pottier J, Pio D, and Mateo R. (2014) ecospat: Spatial ecology miscellaneous methods. R package version 1.

Buisson L, Thuiller W, Casajus N, Lek S, and Grenouillet G. (2010) Uncertainty in ensemble forecasting of species distribution. Global change biology 16:1145-1157. Carbone C, Christie S, Conforti K, Coulson T, Franklin N, Ginsberg J, Griffiths M, Holden J, Kawanishi K, and Kinnaird M. (2001) The use of photographic rates to estimate densities of tigers and other cryptic mammals. Animal conservation 4:75-79. Chaiyarat R, Laohajinda W, Kutintara U, and Nabhitabhata J. (1999) Ecology of the goral (Naemorhedus goral) in Omkoi Wildlife Sanctuary Thailand. Nat Hist Bull Siam Soc 47:191-205. Chang J, Wang Y, Istanbulluoglu E, Bai T, Huang Q, Yang D, and Huang S. (2015) Impact of climate change and human activities on runoff in the Weihe River Basin, China. Quaternary International 380:169-179. Cianfrani C, Le Lay G, Maiorano L, Satizábal HF, Loy A, and Guisan A. (2011) Adapting global conservation strategies to climate change at the European scale: the otter as a flagship species. Biological Conservation 144:2068-2080. Cruz-Cárdenas G, López-Mata L, Villaseñor JL, and Ortiz E. (2014) Potential species distribution modeling and the use of principal component analysis as predictor variables. Revista Mexicana de Biodiversidad 85:189-199. Domíguez-Vega H, Monroy-Vilchis O, Balderas-Valdivia CJ, Gienger C, and Ariano-Sánchez D. 2012 Predicting the potential distribution of the beaded lizard and identification of priority areas for conservation. Journal for Nature Conservation 20:247-253. Drake JM, and Beier JC. (2014) Ecological niche and potential distribution of Anopheles arabiensis in Africa in 2050. Malaria journal 13:213. Dussault C, Courtois R, and Ouellet J-P. (2006) A habitat suitability index model to assess moose habitat selection at multiple spatial scales. Canadian Journal of Forest Research 36:1097-1107. Elith J, Graham CH, Anderson RP, Dudík M, Ferrier S, Guisan A, Hijmans RJ, Huettmann F, Leathwick JR, and Lehmann A. (2006) Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29:129-151. Elith J, and Leathwick JR. (2009a) Species distribution models: ecological explanation and prediction across space and time. Annual review of ecology, evolution, and systematics 40:677-697. Elith J, and Leathwick JR. (2009b) Species distribution models: ecological explanation and prediction across space and time. Annual review of ecology, evolution, and systematics 40. Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, and Yates CJ. (2011) A statistical explanation of MaxEnt for ecologists. Diversity and distributions 17:43-57. Franzmann AW. (1981) Alces alces. Mammalian Species:1-7. Guisan A, and Thuiller W. 2005. Predicting species distribution: offering more than simple habitat models. Ecology letters 8:993-1009. Guisan A, and Zimmermann NE. (2000) Predictive habitat distribution models in ecology. Ecological modelling 135:147-186.

Guo C, Lek S, Ye S, Li W, Liu J, and Li Z. (2015) Uncertainty in ensemble modelling of largescale species distribution: effects from species characteristics and model techniques. Ecological modelling 306:67-75. Hansen MC, Potapov PV, Moore R, Hancher M, Turubanova S, Tyukavina A, Thau D, Stehman S, Goetz S, and Loveland TR. (2013) High-resolution global maps of 21st-century forest cover change. Science 342:850-853. Hijmans RJ, and Elith J. (2013) Species distribution modeling with R. R package version 08-11. Huntley B, Berry PM, Cramer W, and McDonald AP. (1995) Special paper: modelling present and potential future ranges of some European higher plants using climate response surfaces. Journal of Biogeography:967-1001. Iannella M, Cerasoli F, D’Alessandro P, Console G, and Biondi M. (2018) Coupling GIS spatial analysis and Ensemble Niche Modelling to investigate climate change-related threats to the Sicilian pond turtle Emys trinacris, an endangered species from the Mediterranean. PeerJ 6:e4969. Jacobose J. (1993) Climatic records from Northern Areas of Pakistan. Culture Area Karakorum Newsletter 3:13-17. Kabir M, Hameed S, Ali H, Bosso L, Din JU, Bischof R, Redpath S, and Nawaz MA. (2017) Habitat suitability and movement corridors of grey wolf (Canis lupus) in Northern Pakistan. PloS one 12:e0187027. Kindt R. (2018) Ensemble species distribution modelling with transformed suitability values. Environmental Modelling & Software 100:136-145. Li S, McShea WJ, Wang D, Shao L, and Shi X. (2010) The use of infrared‐triggered cameras for surveying phasianids in Sichuan Province, China. Ibis 152:299-309. Liu X, Wu P, Songer M, Cai Q, He X, Zhu Y, and Shao X. (2013) Monitoring wildlife abundance and diversity with infra-red camera traps in Guanyinshan Nature Reserve of Shaanxi Province, China. Ecological Indicators 33:121-128. Luedeling E, Kindt R, Huth NI, and Koenig K. (2014) Agroforestry systems in a changing climate—challenges in projecting future performance. Current Opinion in Environmental Sustainability 6:1-7. Levinsky, I., F. Skov, J. C. Svenning, and C. Rahbek. 2007. Potential impacts of climate change on the distributions and diversity patterns of European mammals. Biodiversity and Conservation 16:3803–3816 Barbet-Massin, M., B. Walther, W. Thuiller, C. Rahbek, and F.Jiguet. 2009. Potential impacts of climate change on the winter distribution of Afro-Palaearctic migrant passerines. Biology Letters 5:248 Luo Z, Jiang Z, and Tang S. (2015) Impacts of climate change on distributions and diversity of ungulates on the Tibetan Plateau. Ecological Applications 25:24-38. Maiorano L, Cheddadi R, Zimmermann N, Pellissier L, Petitpierre B, Pottier J, Laborde H, Hurdu B, Pearman P, and Psomas A. (2013) Building the niche through time: using

13,000 years of data to predict the effects of climate change on three tree species in Europe. Global ecology and biogeography 22:302-317. Marks RB. (2011) China: Its environment and history: Rowman & Littlefield Publishers. Marmion M, Parviainen M, Luoto M, Heikkinen RK, and Thuiller W. (2009) Evaluation of consensus methods in predictive species distribution modelling. Diversity and distributions 15:59-69. Meehl GA, Covey C, Delworth T, Latif M, McAvaney B, Mitchell JF, Stouffer RJ, and Taylor KE. (2007) The WCRP CMIP3 multimodel dataset: A new era in climate change research. Bulletin of the American Meteorological Society 88:1383-1394. Miller J. (2010) Species distribution modeling. Geography Compass 4:490-509. Muñoz-Mendoza C, D'Elía G, Panzera A, Villalobos-Leiva A, Sites JW, and Victoriano PF. (2017) Geography and past climate changes have shaped the evolution of a widespread lizard from the Chilean hotspot. Molecular phylogenetics and evolution 116:157-171. Ogawa-Onishi, Y., P. M. Berry, and N. Tanaka. 2010. Assessing the potential impacts of climate change and their conservation implications in Japan: a case study of conifers. Biological Conservation 143:1728–1736 Parmesan C. (2006) Ecological and evolutionary responses to recent climate change. Annu Rev Ecol Evol Syst 37:637-669. Pearson RG, and Dawson TP. (2003) Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Global ecology and biogeography 12:361-371. Piao R. (1995) Population size and distribution of moose in China. Acta Theriol Sin 15:11-16. Ponder WF, Carter G, Flemons P, and Chapman R. (2001) Evaluation of museum collection data for use in biodiversity assessment. Conservation Biology 15:648-657. Pyke CR, Thomas R, Porter RD, Hellmann JJ, Dukes JS, Lodge DM, and Chavarria G. (2008) Current practices and future opportunities for policy on climate change and invasive species. Conservation Biology 22:585-592. Quintero I, and Wiens JJ. (2013) Rates of projected climate change dramatically exceed past rates of climatic niche evolution among vertebrate species. Ecology letters 16:1095-1103. Renner IW, and Warton DI. (2013) Equivalence of MAXENT and Poisson point process models for species distribution modeling in ecology. Biometrics 69:274-281. Roberts TJ, and Bernhard. (1977) The mammals of Pakistan. Roy P, Ravan SA, Rajadnya N, Das K, Jain A, and Singh S. 1995. Habitat suitability analysis of Nemorhaedus goral–a remote sensing and geographic information system approach. Current Science:685-691. Ruete A, and Leynaud GC. (2015) Goal-oriented evaluation of species distribution models’ accuracy and precision: True Skill Statistic profile and uncertainty maps. PeerJ PrePrints. Sala OE, Chapin FS, Armesto JJ, Berlow E, Bloomfield J, Dirzo R, Huber-Sanwald E, Huenneke LF, Jackson RB, and Kinzig A. (2000) Global biodiversity scenarios for the year 2100. Science 287:1770-1774.

Suryawanshi, K.R., Bhatnagar, Y.V., Mishra, C., 2012. Standardizing the double-observer survey

method

for

estimating

mountain

ungulate

prey

of

the

endangered

snow

leopard. Oecologia, 169(3), 581-590. Sedgeley JA.(2001) Quality of cavity microclimate as a factor influencing selection of maternity roosts by a tree‐dwelling bat, Chalinolobus tuberculatus, in New Zealand. Journal of applied ecology 38:425-438. Shackleton DM. (1997) Wild sheep and goats and their relatives: IUCN. Tateishi R, Hoan NT, Kobayashi T, Alsaaideh B, Tana G, and Phong DX. (2014) Production of global land cover data–GLCNMO2008. Journal of Geography and Geology 6:99. Thomas CD, Cameron A, Green RE, Bakkenes M, Beaumont LJ, Collingham YC, Erasmus BF, De Siqueira MF, Grainger A, and Hannah L. (2004) Extinction risk from climate change. Nature 427:145. Thuiller W. (2004) Patterns and uncertainties of species' range shifts under climate change. Global change biology 10:2020-2027. Tumursukh, L., Suryawanshi, K.R., Mishra, C., McCarthy, T.M., Boldgiv, B., 2015 Status of the mountain ungulate prey of the Endangered snow leopard Panthera uncia in the Tost Local Protected Area, South Gobi, Mongolia. Oryx, 1-6. Thuiller W, Lafourcade B, Engler R, and Araújo MB. (2009) BIOMOD–a platform for ensemble forecasting of species distributions. Ecography 32:369-373. Thuiller, W., O. Broennimann, G. Hughes, J. R. M. Alkemade, G. F. Midgley, and F. Corsi. 2006. Vulnerability of African mammals to anthropogenic climate change under conservative land transformation assumptions. Global Change Biology 12:424–440 Van Vuuren DP, Edmonds J, Kainuma M, Riahi K, Thomson A, Hibbard K, Hurtt GC, Kram T, Krey V, and Lamarque J-F. (2011) The representative concentration pathways: an overview. Climatic change 109:5. Virk A, Sheikh K, and Marwat A. (2003) Northern areas strategy for sustainable development Background Paper: Biodiversity. The World Conservation Union (IUCN), Islamabad. Walther G-R, Post E, Convey P, Menzel A, Parmesan C, Beebee TJ, Fromentin J-M, HoeghGuldberg O, and Bairlein F. (2002) Ecological responses to recent climate change. Nature 416:389. Wang T, Hamann A, Spittlehouse DL, and Murdock TQ. (2012) ClimateWNA—high-resolution spatial climate data for western North America. Journal of Applied Meteorology and Climatology 51:16-29. Warren DL, Glor RE, and Turelli M. (2008) Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62:2868-2883. Wattles DW, and DeStefano S. (2013) Moose habitat in Massachusetts: assessing use at the southern edge of the range. Alces: A Journal Devoted to the Biology and Management of Moose 49:133-147.

Wei J, Zhao Q, Zhao W, and Zhang H. (2018) Predicting the potential distributions of the invasive cycad scale Aulacaspis yasumatsui (Hemiptera: Diaspididae) under different climate change scenarios and the implications for management. PeerJ 6:e4832. Wisz MS, Hijmans R, Li J, Peterson AT, Graham C, Guisan A, and Group NPSDW. (2008) Effects of sample size on the performance of species distribution models. Diversity and distributions 14:763-773. Wogan GO. (2016) Life history traits and niche instability impact accuracy and temporal transferability for historically calibrated distribution models of North American birds. PloS one 11:e0151024. Xu J, Sha T, LI YC, ZHAO ZW, and Yang ZL. (2008) Recombination and genetic differentiation among natural populations of the ectomycorrhizal mushroom Tricholoma matsutake from southwestern China. Molecular Ecology 17:1238-1247. Yang L, Huang M, Zhang R, Lv J, Ren Y, Jiang Z, Zhang W, and Luan X. (2016) Reconstructing the historical distribution of the Amur Leopard (Panthera pardus orientalis) in Northeast China based on historical records. ZooKeys:143. Yates CJ, McNeill A, Elith J, and Midgley GF. (2010) Assessing the impacts of climate change and land transformation on Banksia in the South West Australian Floristic Region. Diversity and distributions 16:187-201. Zhang R, Yang L, Laguardia A, Jiang Z, Huang M, Lv J, Ren Y, Zhang W, and Luan X. (2016) Historical distribution of the otter (Lutra lutra) in north‐east China according to historical records (1950–2014). Aquatic Conservation: Marine and Freshwater Ecosystems 26:602-606. Choudhury, M.R., Deb, P., Singha, H., Chakdar, B., Medhi, M., (2016). Predicting the probable distribution and threat of invasive Mimosa diplotricha Suavalle and Mikania micrantha Kunth in a protected tropical grassland. Ecol. Eng. 97, 23-31. Silvertown, J., (2004). Plant coexistence and the niche. Trends Ecol . Evol. 19, 605611.Smeraldo, S., Di Febbraro, M., Bosso, L., Flaquer, C., Guixé, D., Lisón, F., Meschede, A, Juste, J., Prüger, J., Puig-Montserrat, X., Russo, D., (2018). Ignoring seasonal changes in the ecological niche of non-migratory species may lead to biases in potential distribution models: lessons from bats. Biodiv. Conserv. 27, 2425-2441. Woodward, F.I., (1987). Climate and Plant Distribution. Cambridge University Press, CambridgeAbolmaali, S.M.R., Tarkesh, M., Bashari, H., (2018). MaxEnt modeling for predicting suitable habitats and identifying the effects of climate change on a threatened species, Daphne mucronata, in central Iran. Ecol. Inform. 43, 116-123. Parmesan, C., (2006). Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol. Evol. Syst. 37, 637-669.

Loarie, S.R., Carter, B.E., Hayhoe, K., McMahon, S., Moe, R., Knight, C.A., Ackerly, D.D., (2008). Climate change and the future of California's endemic fl ora. PloS one 3, e2502. Cuena-Lombraña, A., Fois, M., Fenu, G., Cogoni, D., Bacchetta, G., (2018). The impact of climatic variations on the reproductive success of Gentiana lutea L. in a Mediterranean mountain area. Int. J. Biometeorol. 62, 1283-1295.

Table 1. Variables used in species distribution modeling. ID 1 2 3 4 5 6

Describe Annual heat: moisture index Annual Mean Temperature Max Temperature of Warmest Month Annual Precipitation Precipitation of Driest Month Water distance

Name AHM BIO1 BIO5 BIO12 BIO14 Waterdis

Reference Baseline Future WorldClim CCSM4 from WorldClim version1.4 version1.4 WorldClim CCSM4 from WorldClim version1.4 version1.4 WorldClim CCSM4 from WorldClim version1.4 version1.4 WorldClim CCSM4 from WorldClim version1.4 version1.4 WorldClim CCSM4 from WorldClim version1.4 version1.4 Inland water from DIVA-GIS (http://www.divagis.org/gdata)

Figure 1 Species habitat and occurrence point

Figure 2 Model selection by ROC and TSS

Current 40000 goral2070rcp85

goral2030rcp26

30000 20000

goral2050rcp85

goral2050rcp26

10000 0

goral2030rcp85

goral2070rcp26

goral2070rcp45

goral2030rcp45 goral2050rcp45

range change

always suitable

Figure 4: Predicted Habitat and Range for HGG for 2030, 2050 and 2070

Figure 5 Potential distribution of the Goral under climate change.

Figure 3 Species distribution in each time frame in each Rcps

Using an ensemble modelling approach to predict the potential distribution of Himalayan gray goral (Naemorhedus goral bedfordi) in Pakistan Shahid Ahmad1✝, Li Yang1✝, Tauheed Ullah Khan1, Kunyuan Wanghe1 Miaomiao Li1, Xiaofeng Luan1a,* 1 School of Nature Conservation, Beijing Forestry University, NO.35 Tsinghua East Road Haidian District, Beijing, 100083, P. R. China ✝ Joint first authors/these authors contributed equally to this work * Correspondent: Xiaofeng Luan School of Nature Conservation, Beijing Forestry University, NO.35 Tsinghua East Road Haidian District, Beijing, 100083, P. R. China. (E-mail: [email protected]; Tel: +86-13910090393; Fax: +0086-10-62336724)

Conflicts of Interest: The authors declare no conflict of interest.