Implications of landscape genetics and connectivity of snow leopard in the Nepalese Himalayas for its conservation

The snow leopard is one of the most endangered large mammals. Its population, already low, is declining, most likely due to the consequences of human activity, including a reduction in the size and number of suitable habitats. With climate change, habitat loss may escalate, because of an upward shift in the tree line and concomitant loss of the alpine zone, where the snow leopard lives. Migration between suitable areas, therefore, is important because a decline in abundance in these areas may result in inbreeding, fragmentation of populations, reduction in genetic variation due to habitat fragmentation, loss of connectivity, bottlenecks or genetic drift. Here we use our data collected in Nepal to determine the areas suitable for snow leopards, by using habitat suitability maps, and describe the genetic structure of the snow leopard within and between these areas. We also determine the influence of landscape features on the genetic structure of its populations and reveal corridors connecting suitable areas. We conclude that it is necessary to protect these natural corridors to maintain the possibility of snow leopards’ migration between suitable areas, which will enable gene flow between the diminishing populations and thus maintain a viable metapopulation of snow leopards.


Results
Spatial structure of the population. Out of the 268 samples of scat, hairs and urine collected in the three areas studied (Fig. 1), 128 were identified as snow leopard (124 scat and 4 hair samples) based on a mitochondrial DNA species identification.
No microsatellite loci were identified in the twenty snow leopard samples, so they were not included in the following analyses. Thus, finally we had 108 microsatellite samples of snow leopard. After data quality filtering, 63 microsatellite genotypes were obtained, corresponding to 22 individuals according to the identity analysis. All loci were free of errors due to large allele drop out and stuttering, but two loci were suspect because of the presence of null alleles with estimated percentages of 20.27% and 14.45%, respectively. These two loci had higher proportion of homozygotes than expected according to the allele frequency.
Bayesian clustering in Structure clearly attributes the samples from Sagarmatha to spatially separated clusters (Figs. 1A and 1C): one in the northern (N-S) and the other in the south-western (SW-S) part. Samples from Upper Manang (UM) and Lower Mustang (LM) are much more mixed (Figs. 1A and 1B). According to the ΔK parameter in Structure, K = 3 is the most probable number of clusters. The conventional Structure plot can be found in Benesova 37 .
Genetic differences between the three areas studied are analysed in Table 1. Observed heterozygosity is moderate and ranges here between 0.54 and 0.61. Based on fixation index (genetic distance), there is less genetic differentiation in LM and UM, moderate genetic differentiation in SW-S while no differentiation in N-S, however overall there is significance genetic differentiation. See Benesova 37 for details.
Characteristics of the range and results of the connectivity analysis. The habitat suitability map is shown in Fig. 2. The suitable areas are indicated by yellow -green colour in this figure.
In the Maxent model, the receiver operating characteristic (ROC) 38,39 results of the AUC value 38,39 was high for the training data (0.974), indicating that the predictions were excellent and potentially useful. A close look at the Maxent's jackknife test of variable importance 38,39 (Fig. 3) shows that: (i) The most important environmental variables defining the distribution of snow leopard within the area studied were altitude (alt), annual mean temperature (annual_mean_temp), annual precipitation (annual_prec), and distance from roads (roads_main_distance). (ii) Some other factors also seem to play a role: land cover (landcover_main), solar radiation (area_solar) and mean diurnal range of temperature (temp_range). (iii) Omitting the distance from roads from the model causes the largest loss of the explanatory power of the model, compared with omitting any of the remaining variables. (iv) Ruggedness index (rugged) by itself is not useful for estimating the distribution of snow leopard.
Spatial distribution of the most suitable habitats is a relatively narrow and compact belt at altitudes between ca 3500 and 4500 m, when 50% is considered as a threshold for the suitability habitat index (Fig. 4). This belt is widest in regions in western and central Nepal (e.g. in Mustang -see Fig. 1) and narrow in eastern Nepal (Fig. 2). The latter is therefore more vulnerable to fragmentation. Patches of the most suitable habitat are typically southern slopes at an altitude of around 4000 m with a relatively cold and dry climate, shrubs, rocks and open grassland. www.nature.com/scientificreports/    www.nature.com/scientificreports/ Landscape genetic approach and data integration. Habitat connectivity was analysed in order to produce resistance distances between particular localities, for which there were genetic samples and a voltage map indicating the corridors that might be used for dispersal or migration (Fig. 1). Voltage maps reveal regional differences in habitat connectivity, in which highly suitable areas for snow leopards in Central and Eastern Nepal ( There are two possible corridors connecting the western areas (LM, UM) with the eastern one (SNP): (C4) and (C5) (Fig. 1C), which in the westward direction merge and pass through three protected areas: Gaurishankar Conservation Area (GCA), Langtang National Park (LNP) and Manaslu Conservation Area (MCA). However, because of high mountain barrier between these protected areas, the corridor goes partially through the Tibetan Qomolangma National Nature Reserve.
For the whole study area, comparing FST using geographic distances revealed a weak effect of IBD (r = 0.22155; Z = 264,714,021.85418; p < 0.01). A weak effect of IBR is also present (r = 0.15785; Z = 455,729.95660; p < 0.001) as was revealed by comparing FST with resistance distances.

Discussion
We revealed that annual mean temperature, altitude, annual precipitation and distance from roads were the four main factors influencing habitat suitability for snow leopard in Nepal. Bai et al. 36 suggested five main factors (driest quarter, ruggedness, altitude, maximum temperature of the warmest month, and annual mean temperature) in the Qomolangma National Nature Reserve, Li et al. 35 reported two main factors: annual mean temperature and ruggedness in Sanjiangyuan National Nature Reserve, and Aryal et al. 34 reported only one factor: annual mean temperature. Thus, all the studies just mentioned stress climatic factors, especially annual mean temperature and altitude as the main determinants of snow leopard presence.
Our results differ from other ones, as regards ruggedness, which did not have an effect on snow leopard presence in our data. However, ruggedness is usually considered as profitable for snow leopard, as it provides shelter and possibilities to escape 7,15 . It is possible that our results are an artefact of the Maxent modelling, because the scale of the ruggedness layer available may not have been sufficiently fine. The preferred altitudinal range of snow leopard is between ca 3500 and 4500 m. This may be because this approximates the area above the treeline and below the line of maximum plant growth (alpine zone), where the snow leopard's main prey grazes.
According to the genetic results, the gene flow between Manang and Mustang is sufficient to prevent the effects of genetic drift in these two regions. The results of the connectivity analysis reveal three possible ways by which www.nature.com/scientificreports/ snow leopards can move between Manang and Mustang: corridors (C1), (C2) and (C3) -see Fig. 1B. Corridor (C1) does not connect the study sites directly, but could facilitate gene flow between them. It is unlikely that corridor (C2) is used, because it involves climbing great heights: Khatung Kang peak (6484 m) and Tilicho Tal lake (4919 m). Dispersing individuals are thought to be less habitat-specific than residents 40 . However, our camera trap studies indicate the opposite: we photographically captured 5-8 snow leopards at each of the sites (both in UM an LM), but never photographed the same individual at both sites (Shrestha et al. in press). Bayesian clustering in Structure showed that some samples in the Lower Mustang (LM) area were genetically close to those in the south-western part of the Sagarmatha (SW-S), which is approximately 300 km apart. The habitat between them contains many snow leopard habitat patches, as shown in Fig. 1A and three snow leopard populations-those in GCA, LNP and MCA. This means that the study sites are potentially connected via a series of "stepping stones" that could facilitate dispersal among local populations. Thus, there appears to be a certain degree of gene flow between these areas, indicating that a fraction of the population disperses over long distances, which is typical of large carnivores 40 . The average daily distance covered by snow leopards was thought to be 1-7 km 7 . However, recent satellite GPS-based data show that snow leopard can cover even 27-40 km per day 15 and that the average home range is 200-500 km 15,41 . Analyses in Circuitscape reveal two access paths to Sagarmatha, one southern and one western path. According to the results of the cluster analysis there should also be a northern path, probably via the Tibetan plateau, as the samples from northern Sagarmatha (SNP) cluster with those from northern Mustang (LM). Lovari et al. 26 suggest that Sagarmatha was recolonized via the Nangpa Pass from the north, a path not considered in the present study, because we only focused on Nepal.
The results of the landscape genetic analysis indicate a certain degree of gene flow between the two areas studied (LM & UM vs. N-S & SW-S). However, the movement of snow leopards between these areas may be difficult, because (i) the corridor connecting these areas is often narrow (to the west from merging of corridors C4 and C5, suitable habitats consist of narrow strips around the peaks of the mountains, as indicated by the IBD and IBR analyses) -see Because of the discrepancy described in the previous paragraph, it is essential to check the functionality of these corridors and identify others to ensure snow leopard's mobility among suitable areas as soon as possible. If this will not be done, survival of its metapopulation may be jeopardized because of impossibility of betweenpopulations movement.
Both the connectivity and genetic analyses indicate that the populations in N-S and SW-S are genetically isolated from each other. A plausible reason for this is that each of these sub-populations was established from a different source population, when snow leopard was recently re-established in Sagarmatha after more than 40 years of absence 26,42 . The gene flow between these sites might have been limited by the barrier effect of two large rivers (Dudh Khosi and Imja) separating them. Also, the only possible corridor between these sub-populations is narrow and frequently used by tourists on the route to Namche Bazaar. There are hundreds of people on the route during the tourist season 43 , which may be a major factor preventing using it by snow leopard 44 . Thus, because of the short time since the re-establishment and limited gene flow, the genetic differences between the N-S and SW-S might have been preserved.
Similarly to other studies, even here the sightings of snow leopards were very rare. This rarity of sightings in the wild suggests that these animals tend to avoid humans 7,44 , in accord with observed effects of human activities on habitat use by other large predators such as grizzly bears, wolves, and tigers [44][45][46][47][48] . Generally, predators avoid large or frequented roads and trails, especially in areas where hunting or harassment is common 44,45,47,49,50 . Thus, humans may be a substantial determinant of snow leopard presence, in accord with results of our study. The spatial distribution of snow leopard and its prey overlaps with human-modified areas in Nepal. Human development and anthropogenic disturbance, like light pollution or human-associated sounds, change large carnivores' behaviour 51,52 and jeopardize their survival.
To summarize: we show that the genetic structure of snow leopard populations is mainly influenced by the proximity of people and trekking routes used by tourists, which pose a barrier to the dispersal of snow leopards. Topography of the terrain also plays a major role, as it determines the occurrence of suitable habitat for snow leopards in Nepal.

Conclusion
The main message of this paper is the delimitation of the areas with suitable habitat for snow leopards in Nepal, which is presented in Fig. 2, and identification of the main corridors connecting these areas (corridors (C1) -(C5) in Fig. 1). We also determined the influence of landscape on the genetic structure of the population. Isolation could be attributed to either geographic barriers or human disturbances. It is necessary to protect these natural corridors, so that snow leopards can migrate between areas of suitable habitats. This is essential for the gene flow between the diminishing or even disappearing populations and for the survival of a viable population of this rare species. However, maintenance of these corridors may also help to preserve habitat connectivity for other species such as wolves, lynx, golden jackal and red foxes detected in camera traps and carnivore DNA sequence 53 Sites LM and UM are in the Annapurna Conservation Area (ACA), sites N-S and SW-S are in the Sagarmatha National Park (SNP). All four sites were described in our previous paper 4 and therefore we do not repeat their description here. In the Thame valley, no snow leopards were detected in either camera traps or scat samples.

Scat sampling and snow leopard species identification.
In total, 268 putative snow leopard samples  56 for the species identification. The camera trap survey, scat sampling and snow leopard species identification methods were designed for studying abundance of snow leopard and its prey and diet of the former-see Shrestha et al. 4 Leopard Sex ID PCR. The sex identification of the verified snow leopard scat samples was done by testing for the presence of the Y chromosome using primers that amplified an intron of the AMELY gene (i.e., a gene only found on the Y chromosome). The size of this fragment (Y chromosome) is about 200bp 24 .
A PCR reaction of the total volume of 7 μl was prepared containing 3.5 μl of Qiagen 2 × Master Mix Buffer, 0.7 μl of 5X Q solution, 0.05 μl of primers (AMELY-F and AMELY-R) 20 μM, 0.7 μl of distilled water, to which 2 μl of extracted undiluted DNA was added. The PCR reactions were done using the following thermocycling conditions: 95 °C for 15 min; followed by 45 cycles of each 94 °C for 15 s, 55 °C for 30 s and 72 °C for 1 min, with a final extension of 72 °C for 10 min. PCR was run in triplicate; a positive and negative control was maintained. The PCR products were run on a 2% agarose gel, stained with ethidium bromide and visualized under ultraviolet light. The PCR was done in triplicate (3 for each sample), three or two out of three male positives were considered positive. One of three male positives was always repeated to re-confirm the identity. If the repeats were 1/3 for each triplicate, they were considered male based on 6 replicates. Only 3 of 3 negatives were considered to be female. If only one out of six replicates revealed Y chromosome, the samples were regarded as of unconfirmed sex.
Microsatellite based Individual ID. A set of 6 microsatellite loci located on 6 different chromosomes of snow leopard was targeted using sets of the following six fluorescent dye tagged primers in two combinations ( Table 2). Six polymorphic microsatellite loci chosen were sufficient to give the value of PID (Probability of Identity) that was adequate to characterize individuals in a population sharing the same genotype 24 .
PCR reaction was done using a volume of 7 μl composed of 3.5 μl of Qiagen 2X Master Mix buffer, 0.7 μl of 5X Q solution and 2 μl of extracted DNA. 0.88 μl of 20 μM primers were used in the first combination, 0.78 μl of 20 μM primers and 0.02 μl of distilled water were used in the second combination. A multiple tube approach 58 was used for the multiple PCR reactions of each sample, in this case three times. PCR product from both combinations was further diluted to 1:50 and processed using an ABI 31 Analyzer. Three replicates of all the samples were processed. www.nature.com/scientificreports/ Genetic analysis and population structure. After processing the 6-microsatellite based PCR product through capillary electrophoresis using an ABI 31 Analyzer, the labelled fragments by size were separated and retrieved in .fsa files. All these .fsa data files were analysed using different software to assign the allele calls, to determine the genetic diversity and to describe the population genetic structure by Bayesian clustering in Structure (see Benešová 37 for details). Allele calling was done using GeneMarker V1.85 (https ://www.softg eneti cs.com), then consensus genotypes were reconstructed and values adjusted using Autobin software (https :// www6.borde auxaq uitai ne.inra.fr/bioge co_eng-/Scien tific -Produ ction /Compu ter-softw are/Autob in). Presence of genotyping errors due to potential artefacts such as a large allele drop out, stuttering bands etc. was estimated using Micro-Checker 59 .
Occurrence of identical genotypes was analysed using the software Cervus 3.0.7 60 . Basic descriptive parameters of the population genetics were computed using GeneAlEx 6.5 61 . Population genetic structure was described using Bayesian clustering and software Structure 2.3.4 62,63 assuming the clusters differ relative to each other in allele frequency. An analysis was made of the most probable number of clusters and the probability of assigning individual samples to these clusters 63 . Models that presume correlated allelic frequencies and incorporate a biogeographical ancestry analysis were used 64 . Five iterations for each K with 10 6 MCMC steps and burn-in set to 10 5 were performed. The possible number of populations was set in the range of 1 to 10 and for each possible number of populations there were 5 runs. Results were visualized using Structure Selector 65 , the most probable number of clusters was determined using parameter ΔK 66  Habitat suitability and connectivity analyses. Models of habitat suitability indicate both the actual and potential occurrence of a focal species 68 . Beside the definition of suitable or stepping stone areas important for dispersal or migration, outputs of habitat suitability models were used for the construction of resistance surfaces. Inputs needed for habitat suitability models are data on the occurrence of the focal species and relevant environmental variables.
Based on species-habitat associations method by following Gavashelishvili and Lukarevskiy 69 for description of the environmental conditions and preparation of the habitat suitability model, habitat variables related to terrain (topography), climate, habitat (land cover) and effect of human activity (distance to nearest roads) were used based on the published papers 7,15,34,44,53 . The field data were obtained from the Shuttle Radar Topography Mission (SRTM) elevation data and these data were also used to calculate the potential annual solar radiation (Megajoule (mj) cm −2 year −1 ), following Gavashelishvili and Lukarevskiy 69 . The ruggedness data was obtained following Sappington et al. 70 and Riordan et al. 21 . The ruggedness index in the output raster can range from 0 (no terrain variation) to 1 (complete terrain variation). All these layers of environmental variables were prepared in ArcMap 10.6.1. All data were converted into a uniform raster with a 100-m spatial resolution in ASCII format. Habitat analysis and potential distribution modelling were run using Maxent 3.4.1 38,39 and 25% of the occurrence data was used to verify the model. Results of the Maxent model were verified by ROC values whereas models having AUC scores > 0.75 were considered potentially useful and indicate excellent predictive ability 38,39 .
Based on pilot studies (not presented here), the set of best-performing environmental variables was selected, which was then used in the final model. These variables are presented in Table 3.
Inverse raster of habitat suitability provided a resistance surface enabling a further analysis of connectivity. Euclidean distances between localities were calculated using a Geographic Distance Matrix Generator 1.2.3 71 . Landscape connectivity and resistance distances in the area modelled were analysed using the circuitscape Table 3. Environmental variables used in the habitat suitability model. Those that were most closely associated with the presence of snow leopard in the region studied are in bold. www.nature.com/scientificreports/ theory 72 . Software Circuitscape 4.0 73 generated resistance distances between all the localities of the genetic samples, which is the most important output for further integration of geographical and genetic methods. At the same time, voltage maps showing the most probable direction in which individuals moved and consequently the directions of gene flow were produced. The interconnection of the patches where the habitat was suitable for snow leopard was obtained from the habitat model. The most suitable areas were then defined as those for which the habitat model predicts suitability index > 0.5 and an area of at least 10 km 2 .
Landscape genetic approach and data integration. To assess the spatial structure of the population, we used the landscape genetic approach 74 . To analyse the role of particular environmental factors in the spatial ecology of snow leopard, isolation by distance (IBD) 75 and isolation by resistance (IBR) 76,77 were used. Genetic distance measured as F ST 76 was obtained from Genepop 4.2 78,79 . The genetic distance matrix was compared with the geographic and resistance distance matrices using the Mantel test 80 in PASSaGE 2.0 software 81 .
Ethical approval. The research permit was approved by Department and National Parks and Wildlife Conservation and National Trust for Nature Conservation of Nepal.

Data availability
Occurrence locations of snow leopards cannot be presented here due to risk of poaching. We can provide the data to individual researchers upon formal request to either of the authors. www.nature.com/scientificreports/