Genetic analyses reveal population structure and recent decline in leopards (Panthera pardus fusca) across the Indian subcontinent

Research permissions and ethical considerations

All required permissions for our field surveys and biological sampling were provided by the Forest Departments of Uttarakhand (Permit no: 90/5-6), Uttar Pradesh (Permit no: 1127/23-2-12(G) and 1891/23-2-12) and Bihar (Permit no: Wildlife-589). Due to non-invasive nature of sampling, no ethical clearance was required for this study.

Sampling

To detect population structure and past population demography it is important to obtain genetic samples from different leopard habitats all across the study area. In this study, we used leopard genetic data generated from non-invasive samples collected across the Indian subcontinent. We conducted extensive field surveys across the Indian part of Terai-Arc landscape (TAL) covering the north-Indian states of Uttarakhand, Uttar Pradesh and Bihar between 2016–2018. This region has already been studied for large carnivore occupancy using traditional camera trapping as well as field surveys (Johnsingh et al., 2004; Harihar, Pandav & Goyal, 2009; Jhala, Qureshi & Gopal, 2015; Chanchani et al., 2016). We foot surveyed all existing trails covering the entire region to collect faecal samples. Number of trails walked in a particular area was decided based on existing knowledge of leopard presence by the local people and frontline staff members of the sampling team. We collected a total of 778 fresh large carnivore faecal samples. These samples were collected from both inside (n = 469) and outside (n = 309) protected areas from different parts of this landscape. In the field, the samples were judged as large carnivores based on several physical characteristics such as scrape marks, tracks, faecal diameter etc. All faecal samples were collected in wax paper and stored individually in sterile zip-lock bags and stored inside dry, dark boxes in the field for a maximum of two weeks period (Biswas et al., 2019). All samples were collected with GPS locations and were transferred to the laboratory and stored in −20 °C freezers until further processing.

In addition to the north Indian samples collected in this study, we used genetic data previously described in Mondol et al. (2015), representing mostly the Western Ghats and central Indian landscape. The data was earlier used in forensic analyses to assign seized leopard samples to their potential geographic origins in India (Mondol et al., 2015). Out of the 173 individual leopards described in the earlier study, we removed data from related individuals and samples with insufficient data (n = 30) and used the remaining 143 samples for analyses in this study. These samples were collected from the states of Kerala (n = 5), Tamil Nadu (n = 4), Karnataka (n = 53), Andhra Pradesh (n = 3), Madhya Pradesh (n = 12), Maharashtra (n = 46), Gujarat (n = 2), Rajasthan (n = 5), Himachal Pradesh (n = 8), Jharkhand (n = 1), West Bengal (n = 2) and Assam (n = 2), respectively. The sample locations are presented in Fig. 1.

DNA extraction, species and individual identification

For all field-collected faecal samples, DNA extraction was performed using protocols described in Biswas et al. (2019). In brief, each frozen faeces was thawed to room temperature and the upper layer was swabbed twice with Phosphate buffer saline (PBS) saturated sterile cotton applicators (HiMedia). The swabs were lysed with 30 µl of Proteinase K (20 mg/ml) and 300 µl of ATL buffer (Qiagen Inc., Hilden, Germany) overnight at 56 °C, followed by Qiagen DNeasy tissue DNA kit extraction protocol. DNA was eluted twice in 100 µl preheated 1X TE buffer. For every set of samples, extraction negatives were included to monitor possible contaminations.

Species identification was performed using leopard-specific multiplex PCR assay with NADH4 and NADH2 region primers described in Mondol et al. (2014) and cytochrome b primers used in Maroju et al. (2016). PCR reactions were done in 10 µl volumes containing 3.5 µl multiplex buffer mix (Qiagen Inc., Hilden, Germany), 4 µM BSA, 0.2 µM primer mix and 3 µl of scat DNA with conditions including initial denaturation (95 °C for 15 min); 40 cycles of denaturation (94 °C for 30 s), annealing (T_a for 30 s) and extension (72 °C for 35 s); followed by a final extension (72 °C for 10 min). Negative controls were included to monitor possible contamination. Leopard faeces were identified by viewing species-specific bands of 130 bp (NADH4) and 190 bp (NADH2) (Mondol et al., 2014) and 277 bp (cytochrome b) (Maroju et al., 2016) in 2% agarose gel.

For individual identification, we used the same panel of 13 microsatellite loci previously used in Mondol et al. (2014) (Table 1). To generate comparable data with the samples used from earlier study by Mondol et al. (2014) we employed stringent laboratory protocols. All PCR amplifications were performed in 10 µl volumes containing 5 µl Qiagen multiplex PCR buffer mix (QIAGEN Inc., Hilden, Germany), 0.2 µM labelled forward primer (Applied Biosystems, Foster City, CA, USA), 0.2 µM unlabelled reverse primer, 4 µM BSA and 3 µl of the faecal DNA extract. The reactions were performed in an ABI thermocycler with conditions including initial denaturation (94 °C for 15 min); 45 cycles of denaturation (94 °C for 30 s), annealing (T_a for 30 s) and extension (72 °C for 30 s); followed by final extension (72 °C for 30 min). Multiple primers were multiplexed to reduce cost and save DNA (Table 1). PCR negatives were incorporated in all reaction setups to monitor possible contamination. The PCR products were analyzed using an automated ABI 3500XL Bioanalyzer with LIZ 500 size standard (Applied Biosystems, Foster City, CA, USA) and alleles were scored with GENEMAPPER version 4.0 (Softgenetics Inc., State Collage, PA, USA). During data generation from field-collected samples we used one reference sample (genotyped for all loci) from the earlier study for genotyping. As the entire new data is generated along with the reference sample and the alleles were scored along with the reference genotypes, the new data (allele scores) were comparable with earlier data for analyses.

To ensure high quality multi-locus genotypes from faecal samples, we followed a modified multiple-tube approach in combination with quality index analyses (Miquel et al., 2006) as described previously for leopards by Mondol et al. (2009a) and Mondol et al. (2014). All faecal samples were amplified and genotyped four independent times for all the loci. Samples producing identical genotypes for at least three independent amplifications (or a quality index of 0.75 or more) for each loci were considered reliable and used for all further analysis, while the rest were discarded.

Analysis

For each locus, we calculated average amplification success as the percent positive PCR (Broquet & Petit, 2004) after four repeats across all samples. We quantified allelic dropout and false allele rates manually as the number of dropouts or false alleles over the total number of amplifications, respectively (Broquet & Petit, 2004), as well as using MICROCHECKER v 2.2.3 (Van Oosterhout et al., 2004). The false allele frequency is calculated for both homozygous and heterozygous genotypes as the ratio of the number of amplifications having one or more false alleles at a particular locus and the total number of amplifications while allele dropout rate (ADO) is calculated as the ratio between the observed number of amplifications having loss of one allele and the number of positive amplifications of heterozygous individuals.

Post data quality assessment we selected only those samples with good quality data for at least nine or more loci (out of 13) for further analyses. We used the identity analysis module implemented in program CERVUS (Kalinowski, Taper & Marshall, 2007) to identify identical genotypes (or recaptures) by comparing data from all samples. All genetic recaptures were removed from the data set. GIMLET (Valiere, 2002) was used to calculate the PID_(sibs) for all the unique individuals. Following this, any allele having less than 10% frequency across all amplified samples were rechecked for allele confirmation. ARLEQUIN (Excoffier, Laval & Schneider, 2005) was used to determine Hardy Weinberg equilibrium and linkage disequilibrium for all the loci. Finally, to avoid the effects of related individuals in all analyses, we used the program GENECLASS 2.0 (Piry et al., 2004) to select out related individuals in our samples.

To determine the genetic structure of leopards across the Indian subcontinent we used a Bayesian clustering approach implemented in program STRUCTURE (Pritchard, Stephes & Donnelly, 2000; Falush, Stephens & Pritchard, 2003). We performed 10 independent analyses for each K values between one and ten, using 450,000 iterations and a burn-in of 50,000 assuming correlated allele frequencies. The optimal value of K was determined using STRUCTURE HARVESTER web version (Earl & VonHoldt, 2012). Further, we used multivariate analyses approach implemented in program Discriminant Analysis of Principal Component (DAPC) (Jombart, Devillard & Balloux, 2010) to identify genetic clusters in our data. This approach transforms the genetic data into principal components, followed by clustering to define group of individuals with a consideration of minimum within group variation and maximum between group variations among the clusters. The analyses were conducted using adegenet package 2.1.1 in R studio 1.1.453 (R Development Core Team, 2014) where optimal number of clusters were determined through the Bayesian Information Criterion (Jombart, Devillard & Balloux, 2010). Subsequent summary statistics were calculated in ARLEQUIN 3.1 (Excoffier, Laval & Schneider, 2005) and indices of overall genetic differentiation (pairwise F_st) were estimated using GenAlEx version 6.5 (Peakall & Smouse, 2012), dividing the leopard populations according to the STRUCTURE results across the Indian subcontinent. The divisions were based on Q-values (estimated proportions of ancestry) calculated in STRUCTURE, where we used Q > 0.75 as threshold for assigning individuals to a particular population (Mora et al., 2010). Additionally, compression of expected heterozygosity (or G_st) (Nei, 1973) between four leopard sub-populations was calculated in GenAlEx version 6.5 (Peakall & Smouse, 2012). Finally, HP-RARE 1.0 (Kalinowski, 2005) was used to estimate private alleles within each subpopulation.

Demography analyses

Demographic analyses were performed with different genetic subpopulations of leopards based on the results from STRUCTURE analyses. We used a number of different approaches to detect past population demography for leopards. The first two qualitative approaches use summary statistics to detect population size changes, whereas the quantitative approach is a likelihood-based Bayesian algorithm. The summary statistic-based methods used were the Ewens, Watterson, Cornuet and Luikart method implemented in program BOTTLENECK (Cornuet & Luikart, 1996), and the Garza-Williamson index or M ratio (Garza & Williamson, 2001) implemented in program ARLEQUIN 3.1. The quantitative Bayesian approach used was implemented in the program MSVAR 1.3 (Storz & Beaumont, 2002).

Estimation of leopard extinction probability

To understand extinction probability across various biogeographic zones of India we analysed patterns and determinants of leopard occurrence as described in Karanth et al. (2009) and Karanth et al. (2010). In this study, we have just divided the earlier information available for leopards for different genetic subpopulations. We applied a grid-based approach to determine current distribution patterns for leopards, where the selection of grids was based on prior information of leopard presence. This involved collating presence-absence information from more than 100 Indian wildlife experts along with historical information of leopard presence involving hunting locations and other taxidermy and museum records. Each grid cell was an average of 2,818 km² in size and we used data from 1,229 grid cells covering 3,463,322 km² area of the Indian subcontinent. This study applied occupancy modelling to examine the influence of ecological and social covariates on patterns of leopard occupancy. We used a maximum likelihood approach for leopard occupancy in PRESENCE. V.2.0 program (Hines, 2006). Covariates likely to influence leopard distribution modelled included presence and extent of protected areas, land cover-land use characteristics, human cultural tolerance and population density. Data for protected areas was retrieved from the World Database on protected areas (http://www.unep-wcmc.org) and topographic maps. Land cover- land use data were derived from Global Land Cover Facility (2000) and further refined based on Roy et al. (2006) and Joshi et al. (2006). A human tolerance index that characterized different Indian states from most to least tolerant was developed based on knowledge about society-culture, law enforcement, hunting patterns and prior field experiences (for details see Karanth et al., 2009; Karanth et al., 2010). Human population density data were derived from LandScan Global Population Data 2000 (http://www.ornl.gov/gist). Based on existing information on species’ ecology we predicted higher occupancy in protected areas, deciduous-grass-scrub land cover types and lower occupancy in less tolerant states and highly populated areas because of direct competition for food and space (Brashares, Arcese & Sam, 2001; Rangarajan, 2001; Parks & Harcourt, 2002; Karanth et al., 2010). We performed pair-wise correlation tests to screen variables for multicollinearity. The occupancy approach accounts for non-detection of species during surveys and inability to survey some sites (see Karanth et al., 2009; Karanth et al., 2010 for additional details). The probability of extinction was calculated as (1- probability of occurrence) (Karanth et al., 2010). We derived leopard extinction probabilities for three separate major landscapes (Western Ghats, central India and north India) as these regions strongly represented our genetic sampling. These extinction probabilities were compared to the genetically derived estimates.

Individual identification of leopards from north Indian landscape

Of the 778 large carnivore faecal samples collected from TAL, we identified 195 faeces to be of leopard origin (25%) using species-specific PCR assays (Mondol et al., 2014; Maroju et al., 2016). In addition, 457 samples were ascertained to be of tiger (59%) and remaining 126 faecal samples did not produce any result (16%) for either of these large felids, possibly due to poor quality DNA. We amplified 13 microsatellite loci panel on these 195 genetically confirmed leopard faecal samples, and after data validation through multiple repeats generated nine or more loci data from 65 faecal DNA. Subsequently, we identified 56 unique leopard individuals from the 65 samples, whereas nine individuals were ascertained as ‘genetic recaptures’. The mean allelic dropout rate for these loci was found to be 0.05, whereas mean false allele rate for all the 13 loci was 0.04, indicating this 13 loci panel has low genotyping error rates. Amplification success ranged between 41–100% from leopard faecal DNA. None of the loci were found to deviate from the Hardy-Weinberg equilibrium and there were no evidence for strong linkage disequilibrium between any pair of loci. Cumulative PID_sibs and PID_unbiased values were found to be 3.91*10⁻⁶ and 2.73*10⁻¹⁶, respectively, indicating a strong statistical support for unambiguous individual identification. Summary statistics for these samples collected across Terai-Arc landscape is provided in Table 1. We identified 26, 21 and nine unique leopard individuals from the states of Uttarakhand, Uttar Pradesh and Bihar, respectively. As the data generated from north India is comparable to the earlier data, we added this 56 unique leopard data to 143 individual genotypes described in Mondol et al. (2014), and overall 199 unique unrelated leopards were used in subsequent population structure, genetic variation and demography analyses.

Leopard population structure and genetic variation across India

Our sampling strategy targeted countrywide leopard populations to assess population structure and genetic variation. From 199 final unique leopard genotypes we removed four samples representing the eastern and northeast India (n = 2 from the states of West Bengal and Assam each, respectively) from further analyses as they represented inadequate sampling from these regions. Bayesian clustering analysis using 13 microsatellite data from the remaining 195 wild leopard individuals showed four distinct genetic subpopulations (K = 4, see Fig. S1), as presented in Fig. 1. The DAPC analyses identified five different clusters (K = 5) using the Bayesian Information Criterion (Fig. S2). However, out of these genetic clusters two of them were overlapping with each other. Overall, both analyses showed the same pattern of four genetic subpopulations. Majority of the samples showed respective group-specific ancestry, with Western Ghats samples representing the first group (henceforth WG, n = 65), the Deccan Plateau-Semi Arid region forming the second (henceforth DP-SA, n = 66), the samples from Shivalik region covering parts of Himalaya and western parts of upper Gangetic plains making the third group (henceforth SR, n = 38), and finally samples from the Terai region covering eastern part of upper and western part of the lower Gangetic Plains samples forming the fourth one (henceforth TR, n = 26), respectively (Fig. 1). However, small number of samples (n = 18) distributed among the four subpopulations showed mixed ancestry. Subsequent analyses revealed that these leopard subpopulations are genetically differentiated (F_st and G_st) at low, but significant levels (Table 2) for all four populations. The F_st value among these populations ranged between 0.028–0.115, whereas the G_st value between 0.023–0.104 (Table 2).

Analyses and with 13 microsatellite loci among the four genetic subpopulations showed a higher mean number of alleles (NA_WG = 11.77 (S.D. 3.85), NA_DP−SA = 10.46 (S.D. 2.71)) and observed heterozygosity (H_oWG = 0.81 (S.D. 0.08), H_oDP−SA = 0.8 (S.D. 0.08)) in Western Ghats and Deccan Plateau-Semi Arid subpopulations, when compared with samples from Shivalik and Terai region subpopulations (NA_SR = 08.46 (S.D. 2.41), NA_TR = 05.00 (S.D. 1.84) and H_oSR = 0.40 (S.D. 0.14), H_oTR = 0.36 (S.D. 0.28), respectively) (see Table 3 for details). However, the allelic size range values were similar in all populations (Table 3). Western Ghats and Deccan Plateau- Semi Arid subpopulations showed higher number of private alleles (2.38 and 0.85, respectively) when compared to Shivalik and Terai subpopulations (0.46 and 0.15, respectively) (Table 3).

Detection of demographic change

We used microsatellite data to investigate signals of demographic changes in each of the four leopard genetic subpopulations across the subcontinent. Both qualitative approaches, the EWCL and the M-ratio methods indicate signatures of population bottleneck. The EWCL approach implemented in the program BOTTLENECK shows 8-10 loci with heterozygote excess depending on the mutation models used, suggesting a loss of rare alleles through population decline for all four subpopulations. Similarly, the M-ratio approach also shows a low ratio between number of alleles (N_A) and the allelic size range in all four subpopulations (M-ratio_WG-0.37 (S.D. 0.09); M-ratio_DP−SA-0.38 (S.D. 0.09); M-ratio_SR-0.33 (S.D. 0.09); M-ratio_TR-0.29 (S.D. 0.15)), indicating signatures of population bottleneck.

In the quantitative MSVAR approach, models with exponential decline scenarios show consistently that the posterior distributions for log (N0) is always lower than log (N1) for all four subpopulations, indicating population decline for leopards across the subcontinent (Table 4 and Fig. 2). Further quantification revealed that the current effective size is varyingly low (12–25%) than the historical effective size, with Western Ghats, Deccan Plateau-Semi Arid, Shivalik and Terai regions losing approximately 75%, 90%, 90% and 88% of their leopard population, respectively (Table 4 and Fig. 2).

Our subsequent analyses also revealed distributions that suggested recent time of declines in all four populations of leopards (Table 4, Fig. 2). The north Indian subpopulations (Shivalik and Terai) and the Deccan Plateau-Semi Arid population showed the most recent decline occurred about 120–125 years before present, respectively. However, the Western Ghats population indicated potential decline around 200 years ago (Table 4 and Fig. 2).

Leopard occurrence and distribution

We examined the factors influencing leopard distribution at a countrywide scale, where the top ranked model incorporating 28 covariates suggested a wide distribution of habitat types (described in Karanth et al., 2009; Karanth et al., 2010). The model also indicated a positive influence of protected areas, higher cultural tolerance of people and negative influence of higher human population densities and (details in Karanth et al., 2009). Areas with cultivated land, barren areas, deciduous forests and rural–urban were strongly associated with higher leopard occurrence. Naive estimated occupancy was 0.52, whereas model estimated probability of occupancy was significantly higher at 0.68, suggesting that leopards are still widely distributed (Fig. 3) in India compared to most other large mammals (as suggested in Karanth et al., 2010). When compared among the overall three major sub-regions (north India (NI), Deccan Plateau-Semi Arid and Western Ghats), we find that average estimated occurrence was lowest in the north India (Psi_NI = 0.63 ± 0.01, Range: 0.05–1.00, 384 cells) compared to Western Ghats (Psi_WG = 0.83 ± 0.02, Range: 0.23–1.00, 90 cells) and Deccan Plateau-Semi Arid (Psi_DP−SA = 0.79 ± 0.005, Range: 0.25–1.00, 818 cells). Overall, average estimated Psi was 0.74 ± 0.006 (1,292 cells).