A pilot study—genetic diversity and population structure of snow leopards of Gilgit-Baltistan, Pakistan, using molecular techniques

Background The Hindu Kush and Karakoram mountain ranges in Pakistan’s northern areas are a natural habitat of the snow leopard (Panthera uncia syn. Uncia uncia) but the ecological studies on this animal are scarce since it is human shy by nature and lives in difficult mountainous tracts. The pilot study is conducted to exploit the genetic diversity and population structure of the snow leopard in this selected natural habitat of the member of the wildcat family in Pakistan. Method About 50 putative scat samples of snow leopard from five localities of Gilgit-Baltistan (Pakistan) along with a control sample of zoo maintained male snow leopard were collected for comparison. Significant quality and quantity of genomic DNA was extracted from scat samples using combined Zhang–phenol–chloroform method and successful amplification of cytochrome c oxidase I gene (190 bp) using mini-barcode primers, seven simple sequence repeats (SSR) markers and Y-linked AMELY gene (200 bp) was done. Results Cytochrome c oxidase I gene sequencing suggested that 33/50 (66%) scat samples were of snow leopard. AMELY primer suggested that out of 33 amplified samples, 21 (63.63%) scats were from male and 12 (36.36%) from female leopards. Through successful amplification of DNA of 25 out of 33 (75.75%) scat samples using SSR markers, a total of 68 alleles on seven SSR loci were identified, showing low heterozygosity, while high gene flow between population. Discussion The low gene flow rate among the population results in low genetic diversity causing decreased diversification. This affects the adaptability to climatic changes, thus ultimately resulting in decreased population size of the species.


INTRODUCTION
Snow leopard or ounce (Panthera uncia) has evolved in the mountainous ranges of Asia (Bonin et al., 2004). Pakistan's northern mountains are also among its natural habitat, specifically the Gilgit-Baltistan (GB) and Chitral (Khyber Pakhtunkhwa). It is considered an indicator of the health of the ecosystem (O'Brien et al., 1985) and an icon of Kharkoo and Astore (Fig. 1). Each sample was packed separately in centrifuge tubes (50 mL) having silica gel to avoid desiccation, labeled and brought to the laboratory and stored at 4 C till further processing. Samples were subsequently processed for DNA extraction without any delay to obtain good quality of gDNA. A male snow leopard scat sample was provided by Khyber Pakhtunkhwa Wildlife Department (Pakistan), which was used as control.

DNA extraction
Genomic DNA was extracted from scat samples, using three protocols; (1) Zhang Method: Sample (1.5 g) was washed with ethanol, purified by phenol and gDNA was isolated by binding buffer and column filter (Bernevig, Hughes & Zhang, 2006).
(2) Phenol-Chloroform Method: Sample was washed by phosphate buffer saline (PBS), purified by phenol-chloroform and the gDNA was precipitated through ethanol or isopropanol (Wang et al., 2003;Pelaia et al., 2004;Glynn et al., 2004). (3) Combined Zhang and Phenyl-Chloroform Method: Sample (0.25 g) was washed with PBS, bile products were removed by starch and purified with phenol-chloroform, while gDNA was isolated by using binding buffer and filter column. Extracted gDNA was visualized on 1% agarose gel (Sambrook, Fritsch & Maniatis, 1989). Amplification Mini-barcode primers pair having sequence TCCACTAATCACAARGATATTGGTAC and GAAAATCATAATGAAGGCATGAGC (Vaglia et al., 2008) was used to amplify the highly conserved mitochondrial gene, cytochrome c oxidase I (COI), region of extracted gDNA using polymerase chain reaction (PCR). For sex identification, AMELY marker was used (Murphy et al., 1999) after being tested on control male scat sample.

Microsatellite analysis
Confirmed snow leopard samples were genotyped by using seven microsatellites (SSR) markers, that is, PUN82, PUN100, PUN124, PUN132, PUN225, PUN229 and PUN327 (Janečka et al., 2008). Amplified products were visualized on 2% agarose gel. Bivariate data was generated on the basis of pattern of bands, and preceded through GenAlEx to export data file in PopGen format to calculate allelic frequency, genetic distance, heterozygosity and homozygosity (Nei, 1972).

DNA extraction and identification
Combined Zhang-Phenol-Chloroform method was efficient, giving significant extraction success rate, average quantity and quality of gDNA (Table 1; Fig. 2). PCR amplification and sequencing of DNA extracted using 190 bp mini-barcode primer of COI gene (Fig. 3) suggested that 33 out of the total 50 (66%) samples belonged to snow leopard.

Sex determination
Polymerase chain reaction amplification of Y-linked AMELY marker (200 bp) was achieved in control male sample and in 21 field collected samples, suggesting a sex ratio of 1.75:1 male to female (21 males and 12 females, Table 2; Fig. 4). The Chi-square test shows that sex ratio in overall samples of snow leopard scats was not significantly different from 1:1 (p > 0.05).

Genetic analysis
The study could amplify 25 samples out of 33 confirmed snow leopard scat samples at annealing temperature of 55 C for 40 cycles using seven SSR markers. It identified 68 alleles on seven loci with band sizes of 90-220 bp; falling in reasonable proximity with those suggested previously (Janečka et al., 2008). Maximum number (14) of alleles was identified for PUN225 locus and the minimum (6) for PUN82 locus (Table 3). Highest number (44) of alleles were recorded in Shagarthang valley, followed by Thally (33), Astore (27), Basho (26) and the lowest (25) in Kharkoo valley with Shannon indices in different populations ranged between 1.17 and 1.58. Allele frequencies of different SSR loci ranged between 0.02 and 0.75 in different populations (Table 4).    Observed heterozygosity (analyzed on PopGen) was lower compared to expected heterozygosity (Fig. 5). Gene flow (Nm), calculated using Wright's F statistics indices for seven microsatellite primers, reflected high level of gene flow (average 1.65: range 1.18-3.87 for different loci) between populations (Table 5). Nei's diversity indices of genetic distances (Table 6) between populations were low (<1).
Dendrogram constructed using UPGMA, showed low genetic differences and high genetic similarities among three populations (Thally, Kharkoo and Astore).   Shagarthang population exhibited one degree separation from the above three populations, while Basho population appeared as a separate clade (Fig. 6). AMOVA suggested no significant variance among populations indicating no populations' differentiation among regions. Variance contributed by variability was 46% among individuals and that within individuals was 54%, while 0% was observed among populations (Fig. 7).

DISCUSSION
In current study DNA was extracted from 50 putative snow leopard scat samples following three protocols (Taberlet, Waits & Luikart, 1999;Frantzen et al., 1998;Murphy et al., 2002;Wasser et al., 1997). Extraction of gDNA from mucosal layer of scat samples using combined Zhang and phenyl-chloroform method (Cobellis et al., 2004;Glynn et al., 2004;Bernevig, Hughes & Zhang, 2006) yielded significant quantity and quality of DNA, which was sufficient for successful amplification of COI gene/mini-barcode, microsatellite and sex determination primers (Thornton et al., 2008;Valledor et al., 2014). The sample from a control male was amplified by AMELY primer, suggesting the validity of the technique in determination of sex using scat samples of snow leopard. AMELY gene (200 bp) is located on Y-chromosome and hence amplified in males only ( Karmacharya et al., 2011;Janečka et al., 2008;Murphy et al., 1999). Out of total 33 confirmed snow leopard scat samples, 21 scats were of male, whereas 12 samples were of female leopards. The possible reason for less number of female leopards could be the ease of their hunting because they usually remain attached with their cubs, so during hunting they cannot escape easily. Also, during post-parturition the vulnerability of female leopards to mortality factors increases because of an increased ecological and energetic cost of parturition and lactation leading to male-biased sex ratio. Gender-biased infanticide could be another cause of less number of females as stated in literature regarding large carnivores (Oli & Rogers, 1996;Whitman et al., 2004;Balme et al., 2013;Bailey, 2005;Balme, Slotow & Hunter, 2009;Swanepoel et al., 2015). We successfully amplified 75.75% of the field collected scat samples using SSR markers. The lack of amplification in other samples could be due to inhibitors (including bilirubins, bile salts and complex carbohydrates) in the extracted DNA templates (Flekna et al., 2007;Monteiro et al., 2001). The effect of such inhibitors was lowered by using a smaller quantity of the extracted DNA template (Flagstad et al., 2002;Rubin et al., 2000), increasing PCR cycles from 35 to 40 (Frantz et al., 2003;Rubin et al., 2000;Levi et al., 2002;Yekta, hung Shih & Bartel, 2004;Nsubuga et al., 2004;Roon, Waits & Kendall, 2005;Bonin et al., 2004) and using 10× solution "S" as an enhancer (Frackman et al., 1998). Enhancer has the ability to dissolve polar and non-polar compounds resulting in facilitated amplification of DNA (Govinda et al., 2011).
In current study, low genetic variation/heterozygosity (0.33-0.5) was observed in snow leopard population of GB compared to the expected value (0.62-0.75). Low heterozygosity has been reported previously for many populations of snow leopards in Central Asia including Southern Mongolia, Central China, Nepal, north-western India, Pakistan, Tajikistan and Kyrgyzstan (Janjua et al., 2019;Karmacharya et al., 2011;Janečka et al., 2008Janečka et al., , 2011Janečka et al., , 2017. The low genetic diversity in these areas caused by smaller isolated size of the population where higher inbreeding results in increased homozygosity and genetic fixation or loss of unexploited genetic potentials (Kotzé & Muller, 1994). It can also be attributed to bottleneck effect, genetic drift and inbreeding depression which sometimes lead to expression of deleterious recessives alleles resulting in lowered survival rate of species (Lindenmayer et al., 1993;Franklin & Frankham, 1998;O'Brien et al., 1985). Small population size resulting in inbreeding depression reduced the chances of adaptation to environmental changes (Barone et al., 1994;Crooks, 2002).

CONCLUSION
The findings of genetic diversity, population structure and function played a vital role in formulation of conservational strategies. The low genetic diversity in snow leopard populations of the study leads to the conclusion that the gene flow among the populations is too low and the genetic diversification of the animal is not enough to aptly adapt to the environmental changes which would not result in the efficient promotion of this species. Therefore, proper planning and management in the protected and non-protected areas is required to get the output. The illegal hunting, poaching and human-animal conflicts have to be stopped by arranging proper protective and safety measures to the natural habitat of the inhabitant species in order to maintain the effective population size of the threatened wild fauna.

Author Contributions
Samreen Aruge performed the experiments, analyzed the data, prepared figures and/or tables. Hafsa Batool performed the experiments, prepared figures and/or tables. Fida M. Khan analyzed the data, contributed reagents/materials/analysis tools, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft, collected samples. Fakhar-i-Abbas conceived and designed the experiments, contributed reagents/ materials/analysis tools, prepared figures and/or tables, approved the final draft, collected samples. Safia Janjua conceived and designed the experiments, approved the final draft.

Data Availability
The following information was supplied regarding data availability: Raw data is available in Table S1. The raw data consists of the allele frequency of SSR markers at different loci in snow leopard population samples analyzed by POP Gen.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.7672#supplemental-information.