Genetic insights to assist management of the Critically Endangered hangul Cervus hanglu hanglu in the Kashmir Himalaya

Abstract

The hangul Cervus hanglu hanglu, a Critically Endangered mountain ungulate of Jammu and Kashmir, India, faces the imminent threat of population loss and extinction. Effective management of its largest viable population in Dachigam National Park in the Kashmir Himalaya requires reliable demographic information. Using 14 microsatellite markers we identified 293 individuals (208 females and 85 males) through faecal analysis, and generated data on the genetic status and population size of the hangul in its winter habitat. The mean expected and observed heterozygosities of 0.62 and 0.59 are comparable to those of several red deer Cervus elaphus populations elsewhere. The effective population sizes were 46.3 and 93.7 when the frequencies of rare alleles were considered to be 0.050 and 0.010, respectively. The average mean kinship of the population was 0.34, and there was no evidence of a recent bottleneck event. In genetic mark–recapture analysis the best model included an effect of sex on both detection and recapture probabilities. Detection of males was highest in November, coinciding with the hangul breeding season, whereas detection of females was highest in December. Our estimate of the hangul population using genetic mark–recapture with bootstrapping was 394 individuals. To our knowledge, this is the first study to use genetic data to estimate the population of the hangul. It will guide future studies of this subspecies and also serve as an impetus for identifying founder animals for captive breeding, and for connecting the population in Dachigam National Park with the other small, isolated populations to ensure the long-term survival of this subspecies.

Introduction

The Anthropocene is characterized by unprecedented levels of biodiversity loss and defaunation (Dirzo et al., 2014). Vertebrate populations have declined by a mean of 60% since 1970, and 96% of all mammalian biomass consists of humans and livestock (Turvey & Crees, 2019). Large-bodied herbivores are more susceptible to extinction not only because of their slow life histories and low fecundity rates, but also because they are disproportionately hunted for meat and other resources, and frequently compete with livestock for habitat (Cardillo et al., 2005; Ripple et al., 2015). Mountain ungulates are at particular risk as they depend on specific habitat conditions and migrate seasonally between elevations in search of suitable pastures and secure landscapes. Small and isolated populations are vulnerable to stochastic environmental events and are prone to local extinctions (Grøtan et al., 2008), and they are also susceptible to genetic drift that reduces genetic diversity, increases inbreeding over time and compromises their ability to persist in changing environments (Poirier et al., 2019).

The hangul Cervus hanglu hanglu was once considered a subspecies of the red deer Cervus elaphus complex along with two other subspecies, the Yarkand–Tarim red deer Cervus elaphus yarkandensis and the Bukhara red deer Cervus elaphus bactrianus of Central Asia (Sun et al., 2008; Lorenzini & Garofalo, 2015; Meiri et al., 2017). The hangul is categorized as Critically Endangered on the IUCN Red List (Brook et al., 2017), and as a Schedule 1 species in the Indian WildLife (Protection) Act, 1972. The hangul was once widely distributed across Kashmir, the Chenab Valley of Jammu, and Himachal Pradesh. However, major range contractions have confined the subspecies to the Kashmir Himalaya (Thakur et al., 2015). Like many other mountain ungulates, the hangul is a seasonal migrant, moving between lower and higher elevations in winter and summer, respectively (Sharma et al., 2010; Ahmad et al., 2016). The hangul feeds mostly on dicotyledonous trees and shrubs, except in spring and summer when forbs and grasses form a major part of its diet (Shah et al., 2009; Sharma et al., 2010; Ahmad et al., 2016). It is a seasonal breeder, with a few observational and anecdotal records suggesting that males rut during September–November, when they guard large female harems, and break away from the groups immediately after the mating season.

There were 3,000–5,000 hangul at the beginning of the 20th century (Qureshi et al., 2009; Thakur et al., 2015). However, according to more recent census reports from the Department of Wildlife Protection, Jammu and Kashmir (2019, 2021), only c. 200 individuals survive in the wild, with a skewed sex ratio of 12.5–15.0 males to 100 females. The only viable population is in Dachigam National Park, along with small populations in Shikargah, Tral Wildlife Sanctuary and Overa Wildlife Sanctuary in Jammu and Kashmir. A few smaller populations are also known from Wangat–Naranag, Chandaji–Diver–Lolab and Overa–Aru (Qureshi et al., 2009; Kaul et al., 2018). The lack of connectivity between these once-contiguous populations has also led to a likely decrease of genetic diversity (Thakur et al., 2013). The decline of the hangul has been driven by anthropogenic factors such as poaching, habitat fragmentation and degradation, competition for space and forage with livestock populations in their summer habitat, predation by guard dogs and disease transmission from livestock (Gee, 1965; Schaller, 1969; Qureshi et al., 2009).

There is thus an urgent need for data to guide effective management of the hangul population in the Kashmir Himalaya so that conservation interventions can restore or maintain genetic diversity and ensure long-term survival of the population. Here we estimate the genetic diversity, population size and sex ratio of the hangul in Dachigam National Park and discuss the implications of our findings for effective conservation of the subspecies.

Study area

We conducted this study in Dachigam National Park, c. 21 km north-east of Srinagar in Jammu and Kashmir, India, during November 2019–February 2020 (Fig. 1). The Kashmir Himalaya falls within the biogeographical unit of the north-western Himalaya, which supports 11 large ungulate species. Dachigam National Park lies in the Zabarwan mountain range, which is part of the Zanskar Himalaya of the north-western branch of the central Himalayan axis. The Park encompasses an area of 144 km², with an elevational range of 1,600–4,400 m. It comprises lower Dachigam (69 km²; 1,600–3,000 m), which is the winter habitat of the hangul, and upper Dachigam (75 km²; 3,000–4,400 m), which is its summer habitat. Lower Dachigam contains several army settlements and guest houses, and the Park boundary is surrounded by a number of villages. Upper Dachigam is visited by livestock herders in summer, during June–September. In addition to the hangul, the large mammals of the area include the common leopard Panthera pardus, Himalayan black bear Ursus thibetanus, red fox Vulpes vulpes, Himalayan yellow-throated marten Martes flavigula and Kashmir musk deer Moschus cupreus.

Fig. 1 The study area and locations of the 11 survey trails in Dachigam National Park, Jammu & Kashmir, India. Trail numbers correspond to those in Table 1. (Readers of the printed journal are referred to the online article for a colour version of this figure.)

Based on temperature and phenological patterns of vegetation in the region, spring is April–May, summer is June–September, autumn is October–November and winter is December–March. The temperature ranges from a minimum of −7 °C in winter to a maximum of 30 °C in summer, with annual precipitation of 628–932 mm. The area is snow-covered from December until March or April. Above 3,000 m snow begins to accumulate in November, with heavy snowfall typically during December–February, and this snow usually melts during April–May.

Methods

Study design

We surveyed the hangul population in its winter habitat (i.e. lower Dachigam), because at this time the population is restricted to a small and accessible area and we would expect the population to have a smaller home range than in spring or summer (as reported for other red deer species; Georgii, 1980; Kleveland, 2007), improving our chances of sampling a large proportion of the population. We marked 11 permanent trails of 1.0–3.6 km, encompassing two habitat types (Table 1, Fig. 1). Trails were at least 2.5 km apart, to maximize our chances of capturing a greater number of individuals. We used these trails to collect faecal pellets and associated relevant information.

Table 1 Details of the 11 permanent survey trails marked in lower Dachigam National Park, Jammu and Kashmir, India (Fig. 1).

Faecal sample collection

We monitored the 11 trails (Fig. 1) once per month during November 2019–February 2020 to collect fresh hangul faecal pellets. We identified hangul faecal pellets visually. There was negligible chance of misidentification as the hangul is the only large herbivore in the area in significant numbers. Musk deer, the only other large ungulate in the region, exists in low numbers, and their pellets are distinctly smaller and thinner, with a tapering end, than those of the hangul. In addition, the musk deer is territorial, and its pellets are mostly found in piles, in contrast to the individual pellet groups of the hangul. We later verified the species identity of pellets through molecular analysis. We collected all fresh faecal samples encountered on a trail, and from each pellet group we collected 3–4 pellets in 90% ethanol in 50 ml vials, which we labelled and stored at −20 °C until analysis. We gave each sample a unique identification number, and recorded its location and elevation, and the date and time of collection.

DNA extraction and genetic profiling

We scraped the outer layers of faecal pellets to obtain DNA of the target species from shed rectal epithelial cells. We extracted DNA using the NucleoSpin DNA stool kit (Macherey-Nagel, Dueren, Germany) according to the manufacturer's instructions, except we employed longer incubation periods to improve DNA yields. We conducted the initial processing of pellets and the DNA extraction in separate facilities, to avoid contamination. We assessed the quality and quantity of DNA from the samples spectroscopically, and we used samples with 10–20 ng of DNA/μl directly, whereas we diluted samples with higher yields to 10 ng/μl for PCR.

We designed species-specific molecular markers to confirm the identity of hangul pellet samples. Unlike previous work (Thakur et al., 2015), we aimed to develop a simple and accurate assay to confirm species identity without the need for an additional sequencing step. We aligned mitochondrial cytochrome b sequences of several ungulate species (the hangul, elk Cervus canadensis, domestic goat Capra hircus, domestic sheep Ovis aries, fallow deer Dama dama, sika deer Cervus nippon, sambar deer Rusa unicolor, musk deer, Himalayan serow Capricornis sumatraensis thar and spotted deer Axis axis) using MEGA X 10.2.4 (Kumar et al., 2018) to identify sequence variations specific to the hangul (Supplementary Fig. 1). Based on these point variations we designed two primers: HSS-F: 5′-CTTCTAGTACTCGTCTTAAC-3′ and HSS-R: 5′-CAGGTTTAATATGAGGGAGC-3′. We confirmed primer specificity with Primer-BLAST (Ye et al., 2012) and conducted in silico PCR in AmplifX 2.0.7 (Jullien, 2020) before analysing all samples.

We screened in triplicate 14 microsatellite loci (CSSM22, CSSM16, CSSM19, BM1818, Cer14, ILSTS06, CSSM66, Haut14, CSPS115, INRA35, ETH225, CSRM60, BM888, BM4208) reported to be polymorphic in European red deer (Kuehn et al., 2003; Zhou et al., 2015). We selected these markers based on their reported polymorphisms, chromosome locations and amplicon sizes in multiple studies on red deer (Bishop et al., 1994; Moore et al., 1994; DeWoody et al., 1995; Thieven et al., 1995; Kühn et al., 1996; Vanessa et al., 2014).

We identified the sex of putative individuals using markers designed to target variable regions in the amelogenin gene on the X and Y chromosomes of red deer (Pfeiffer & Brenig, 2005). See Supplementary Material 1 for all PCR details.

Molecular data analyses

We scored microsatellite results with GeneMapper 5.0 (Applied Biosystems, Waltham, USA). We used CERVUS 3.0.7 for allele frequency and identity analyses to find matching pairs of genotypes and to determine probability of both individual and sibling identity (Kalinowski et al., 2007). We further analysed samples that worked at a minimum of seven loci. We allowed for mismatches in only two loci to rule out misidentification as a result of genotyping errors. We manually re-examined mismatching genotypes to rule out scoring or entry errors. We analysed PCR success rate and accuracy of microsatellite amplification at all 14 loci for 348 putative genotypes. We calculated PCR success rate as the proportion of positive PCRs to the total PCRs set up. We also calculated allelic dropout and false allele frequencies at each locus (Supplementary Table 1). We determined deviations from the Hardy–Weinberg equilibrium with both CERVUS and GENEPOP (Raymond & Rousset, 1995; Rousset, 2008) with default values of Markov chain parameters, and we estimated pairwise linkage disequilibrium using ARLEQUIN 3.5.1.2 (Excoffier & Lischer, 2010). We used NeEstimator 2.01 (Do et al., 2014) to calculate effective population size (N_e) based on the linkage disequilibrium model of random mating, with the critical value or the lowest allele frequency set at 0.050 and 0.010 to screen for rare alleles. Effective population size is the genetic equivalent of a population census, and it represents the number of breeding individuals in an idealized population that contribute genetically to the next generation (Wright, 1931, 1938; Waits & Storfer, 2016). Rare allele screening using critical values helps evaluate the effects of low-frequency alleles on estimates of effective size (Waples, 2006). We calculated the number of observed alleles, observed and expected heterozygosities and polymorphic information content with MolKin 3.0 (Gutiérrez et al., 2005). We also used MolKin to determine mean kinship between individuals. Kinship or co-ancestry between individuals i and j (k_ij) is the probability that two alleles at a locus taken at random, one from each individual, are identical by descent (Wright, 1969; Frankham et al., 2010, 2017). The mean kinship of an individual i (mk_i) is the mean of all kinships between that individual and all living individuals in the population including itself:

$$mk_i = \mathop \sum \limits_{\,j = 1}^{\rm N} k_{ij}/{\rm N},$$

where N is the number of individuals in the population (Frankham et al., 2017). We also assessed the possibility of an historical genetic bottleneck, using BOTTLENECK 1.2.01 (Piry et al., 1999). This software uses allele frequency data to analyse heterozygosity excess at Hardy–Weinberg equilibrium and mutation drift equilibrium. A higher expected heterozygosity value than expected equilibrium heterozygosity (H_e > H_eq) at the majority of the loci indicates a recent bottleneck event. The allele frequency distribution (mode-shift indicator) test, which is qualitative, gives a graphical representation of data to demarcate stable populations and those that have undergone a bottleneck (Luikart et al., 1998).

Population estimation

We used encounter histories of individuals identified through their genetic profiles to estimate the hangul population in all of the months sampled during November 2019–February 2020 in Dachigam National Park. We surveyed all 11 permanent trails once each month over 5 days in November, 6 days in December, 7 days in January and 6 days in February (24 days in total). The gap between trail visits within each month was < 2 days, except in November when one of the gaps between visits exceeded 6 days because of severe snowfall during that time. There was a gap of at least 15 days between trail surveys in consecutive months. We organized the encounter data for analysis as per the robust design framework in Program MARK (White & Burnham, 1999), details of which are provided in Supplementary Material 2.

We estimated abundance using closed population model assumptions. We built nine models and selected the most parsimonious, using the Akaike information criterion corrected for small sample size (AICc; Burnham & Anderson, 2002). We performed model averaging using MARK for those models that had ΔAICc < 2.

We also estimated the population by bootstrapping the dataset. We used the packages rich 1.0.1 (Rossie, 2016) and vegan 2.5-7 (Oksanen et al., 2020) in R 4.1.3 (R Core Team, 2013). We created a matrix of sampling days (rows) and individuals sampled on each day (columns), and we bootstrapped the data 1,000 times with samples drawn at random from the dataset. This procedure considered all of the captured individuals as a data universe for bootstrapping and estimated the population and the standard error of the estimate by randomly subsampling the data universe.

Results

We collected 547 fresh faecal samples during November 2019–February 2020, and extracted a usable quantity of DNA (> 20 ng/μl) from all samples. We confirmed species identification of these samples with our new primer set, HSS-F and HSS-R, which yielded a hangul-specific amplicon of 112 bp (Supplementary Fig. 2). Of these, 348 samples could be genotyped at seven or more microsatellite loci. The probability of identity of unrelated individuals was 7.7 × 10⁻¹², whereas the sibling probability of identity was 3.4 × 10⁻⁵, with 14 loci. Null allele frequencies ranged from −0.25 to 0.16, indicating a high proportion of heterozygotes in the data. PCR success ranged from 84% (locus INRA35) to 96% (Cer14). INRA35 also showed the highest rate of false alleles (18%) and allelic dropouts (32%). However, the data on loci CSSM22, BM1818 and Cer14 did not contain false alleles, and Haut14 had the lowest dropout rate (4%). We used the dropout rate at each locus to calculate the number of PCR repetitions needed to obtain results with > 99% certainty in the homozygous condition. This level of certainty could be achieved with three repeats for most of the loci except for ILSTS06, INRA35 and ETH225, which required four repeats (Supplementary Table 1).

We resolved 293 individuals using identity analysis in CERVUS 3.0.7, of which 43 individuals were recaptured at least once in our dataset (Table 2). Following this analysis, we attempted to fill any gaps in the final 293 genotypes, and overall they were 71% complete. Of these 293 individuals, 208 were female and 85 were male. Allele frequency analysis of these genotypes yielded mean expected and observed heterozygosities of 0.62 and 0.59, respectively (Table 3). The number of alleles ranged from 2 (locus CSSM22) to 12 (BM888), with a mean of 6.29 alleles. Following Bonferroni correction, most of the loci did not show significant deviation from Hardy–Weinberg equilibrium, except for Cer14 and Haut14. Seven pairs of microsatellite loci of 91 possible combinations showed highly significant levels of association. The effective population size (N_e), when the rare allele frequency in the data was assumed to be 0.05, was 46.3 (34.7–62.8), and when the frequency was assumed to be 0.01, this value was 93.7 (72.2−125.8). Mean kinship in the population ranged from 0.21 to 0.50, with a mean overall value of 0.34. We obtained a normal L-shaped distribution in the mode-shift indicator test, which rules out any recent bottleneck events in the population (Supplementary Fig. 3).

Table 2 Details of hangul Cervus hanglu hanglu faecal samples collected in Dachigam National Park during November 2019–February 2020, with the number of genotypes and genotyped individuals, and the number of these recaptured.

Table 3 Statistical indices of genetic diversity of the hangul population in Dachigam National Park obtained from 14 polymorphic microsatellite loci.

In the mark–recapture analysis, the best model included an effect of sex on both detection (p) and recapture probability (c), and p was not equal to c in all of the months for both females and males (Supplementary Table 2). The detection probability of male hanguls was highest in November, during the hangul breeding season (Table 4). The mean detection probabilities p of both male and female hanguls in December were 0.08 (95% CI 0.005–0.061) and 0.25 (95% CI 0.16–0.38), respectively. In January, however, there was a steep decline in detection probabilities for both males and females, which only slightly improved in February (Supplementary Fig. 4). Recapture probability was close to zero in all months for both sexes (Table 4). This could be because of the sampling design, in which we avoided collecting old faecal samples from transects. Although we were not able to estimate the population size with this analysis, estimates in November for males, and in December and February for females are similar to the number of males and females genotyped in the respective months (Table 2, Supplementary Table 3).

Table 4 Model-averaged estimates of detection and recapture probabilities for male and female hanguls across the four survey months in Dachigam National Park.

The bootstrapping method gave a population estimate of 394 ± SE 26 individuals (95% CI 293–445). Other estimators, as in the genetic mark–recapture analysis above, would probably be highly biased as the sample pool remains unsaturated. The bootstrapping method provides a better approximation than mark–recapture analysis and can be considered a conservative estimate of the population in the sampled area.

Discussion

There are several caveats to our sample collection. We surveyed in winter to improve our chances of capturing a large proportion of the hangul population. However, the harsh winter months in the mountains constrain data collection. We could not monitor our trails at the scheduled time during spells of heavy snowfall, and there were instances when some trails could not be accessed at all within a month. Fresh faecal samples were often buried under snow, thereby reducing our encounter rate. During heavy snowfall hanguls might have moved to elevations below our permanent trails, further affecting faecal collection. Nevertheless, our data collection and DNA analysis is a baseline for planning future studies.

The genetic diversity of the Critically Endangered hangul population (H_e = 0.62) in Dachigam National Park is comparable to that of several red deer populations (categorized as Least Concern on the IUCN Red List) across Europe and parts of Asia (Kuehn et al., 2003; Dellicour et al., 2011; Frantz et al., 2012; Šprem et al., 2013; Hoffmann et al., 2016; Zachos et al., 2016; Galarza et al., 2017; Edelhoff et al., 2020), with overlapping microsatellite markers in many of the studies. This finding provides hope for the survival and expansion of the largest extant hangul population in this region. Furthermore, the levels of genetic diversity explain why signatures of a population bottleneck, inbreeding or a historical population crash were not detected in earlier studies (references in Ahmad et al., 2023). However, the linkage disequilibrium results and Hardy–Weinberg equilibrium analyses in our study indicate possible early signs of genetic drift, population sub-structuring and even isolation (Johansson et al., 2005; Li & Merilä, 2010; Kiselyova et al., 2014). The average mean kinship for a population is the expected inbreeding coefficient in the next generation with random mating. It is also directly related to the proportion of genetic diversity lost since time 0. Consequently, if kinship is minimized, inbreeding in the next generation is minimized and proportional heterozygosity retained is maximized (Frankham et al., 2017). Mean kinship in the hangul population in Dachigam National Park is 0.21–0.50, and the population average of this mean kinship is 0.34. This could be a valuable parameter for identifying individuals or populations with low kinship values when planning captive breeding or translocation programmes (Frankham et al., 2017). Average mean kinship is considered more informative than F-statistic values such as F_ST (fixation index) and F_IS (inbreeding coefficient) in understanding the level of inbreeding within and between populations. Effective population size (N_e) is generally smaller than census population size (Frankham, 1995), and in this case N_e is 46.3–93.7 for the minimum of 293 identified individuals. Factors contributing to smaller N_e include an unequal sex ratio, assortative mating and genetic drift, all of which can lead to declines in genetic diversity in small, isolated populations.

Of the 293 hangul individuals identified, 208 were female and 85 male, a sex ratio of c. 2.5 females to 1 male. However, earlier census reports showed highly skewed adult sex ratios (Department of Wildlife Protection, Jammu and Kashmir, 2019, 2021). There are several possible reasons for these differing results. We collected faecal samples towards the end of the rutting season (Bhat et al., 2009) and through peak winter when adult males are expected to be living with or in close proximity to female herds. However, the censuses are usually carried out in the spring and early summer months (Qureshi et al., 2009; Department of Wildlife Protection, Jammu and Kashmir, 2019, 2021), by which time adult males have shed their antlers (J & K ENVIS Newsletter, 2016), separated from herds and dispersed over wide areas. These factors may have confounded accurate identification of sex through direct sightings or photography. However, the sex ratio in our present study should be interpreted cautiously, as we could not distinguish young from adults through DNA analysis of faecal samples.

This is the first faecal DNA-based genetic mark–recapture study of the hangul. Although there are no accurate data on former numbers of hangul, populations of 3,000–4,000 in the early 1900s suffered severe declines from the 1960s onwards (Gee, 1965; Schaller 1969; Holloway 1970; Kaul et al., 2018; Srivastava & Vasudevan, 2021). In non-invasive mark–recapture studies, trap shyness or attraction to traps is generally minimal (Boulanger et al., 2006). However, the terrain we sampled received variable snowfall across the months of the study. After heavy snowfall it is difficult to find faecal pellets, and this was reflected in the low capture and recapture probabilities we observed in January for both females and males. However, an effect of sex was also apparent in the estimates for November, in which the male detection probability was much higher than that of females. The hangul rutting season starts in late September and runs until mid November (Bhat et al., 2009), and thus our detection of males may have been higher than recorded previously because our study took place during this period. Detection of females was highest in December, indicating the suitability of this period for conducting surveys. Eventually males disperse from the herds, which could be why there were lower detections of males in February. Consistent trap shyness (i.e. the difficulty of capturing and recapturing individuals) could also be because we collected only fresh faecal pellets for genotyping: potential captures or recaptures might have been missed if hanguls passed across the trails more than a few days before the survey began.

Our findings provide a basis for developing robust and accurate methods for estimating populations of mountain ungulates and the hangul in particular. November and December appear to be the best months in which to monitor the hangul population in lower Dachigam. Survey design for genetic mark–recapture can be improved to estimate population size reliably and can also be done in combination with advanced camera-trapping methods (Pal et al., 2021). It is also important to use genetic studies to evaluate the degree of connectivity of the Dachigam population with the other small, isolated hangul populations in the Kashmir Himalaya and to take steps to improve gene flow and revive these populations. In addition, a genetically diverse founder population needs to be established in captivity for future restocking and strengthening of the wild populations.