Evolutionary insights on critically endangered Kashmir red deer or hangul (Cervus hanglu hanglu) through a mitogenomic lens

Research permits and ethical considerations

All required research permissions towards fieldwork and biological sampling were accorded by the Chief Wildlife Warden, Department of Wildlife Protection, Jammu and Kashmir (Permit Nos. WLP/TECH/116-19, WLP/TECH/871-74/2013, and WLP/F-101/Res/2019-20/424-27). Blood samples (n = 2) were collected as part of an ongoing study titled “Long-term conservation plan for Hangul: movement pattern study using satellite telemetry”, where invasive sampling was conducted during radio-collaring operations at Dachigam National Park, Kashmir in 2013 and 2020, respectively. Ethical permissions were provided by the Sher-e-Kashmir University of Agricultural Sciences and Technology of Kashmir, Srinagar, Jammu and Kashmir, India (letter no. AU/FVS/PS-23/22/8553).

Sampling

Apparently healthy adult Hangul (n = 2, 1 M:1 F) individuals captured for satellite collaring were sampled in this study. The male was captured in 2013 (estimated weight of ~150 kg) whereas the female was captured in 2020 (estimated weight of ~120 kg). In both cases, the animals were lured into an open forest patch using green fodder and salt and were chemically immobilized using a combination of medetomidine (10 mg/ml, Wildlife Pharmaceuticals Inc., Fort Collins, CO, USA) and ketamine hydrochloride (100 mg/ml, Troy laboratories Pvt. Ltd., Glendenning, NSW, Australia) at dose rates of 90 μg kg⁻¹ and 1.45 mg kg⁻¹ for male and 120 μg kg⁻¹ and 2.08 mg kg⁻¹ for female, respectively (Jalanka, 1988; Jalanka & Roeken, 1990). Drug delivery was carried out remotely from a hide using Dan-inject syringe projector (Model-IM) from a distance of 30–40 m. The animals were approached after achieving sedation safe for handling, and blood sampling was conducted. Blood was collected aseptically (5 ml from each animal) by jugular vein puncture in EDTA BD Vacutainer (BD Diagnostics, Franklin Lakes, NJ, USA) and transported to the laboratory within 2 h and stored at −20 °C till further processing. Drug reversal was assisted by administering 5 mg/ml atipamezole (Dechra, Northwich, United Kingdom) intramuscularly at a dose rate of 0.2 mg kg⁻¹ body weight following completion of necessary procedures.

DNA extraction and mitogenome sequencing

Mitogenome sequencing of Hangul was performed using next-generation sequencing (NGS) approach. Genomic DNA was extracted from two blood samples using MN NucleoSpin Blood kit (Macherey-Nagel, Duren, Germany) following manufacturers’ instructions. DNA quality was checked through 1% agarose gel electrophoresis and Qubit 3.0 fluorometer (Invitrogen, Life Technologies, Waltham, MA, USA). Following this, paired-end libraries were prepared using NEBNext Ultra DNA Library Prep kit (New England Biolabs, Ipswich, MA, USA) where the genomic DNA was enzymatically fragmented, end-polished, A-tailed, and ligated with full length Illumina sequencing adapters. The PCR amplification was conducted using universal P5 and indexed P7 oligos with following conditions: initial denaturation at 98 °C for 30 s, followed by four cycles of denaturation (98 °C for 10 s), annealing (65 °C for 75 s) and final extension (65 °C for 5 min). The constructed libraries were checked for purity, size, and concentration using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Finally, the qualified libraries (2 × 150 bp) were sequenced on Illumina NovaSeq6000 (Illumina, San Diego, CA, USA) platform at the Next Generation Genomics Facility of Nucleome Informatics Pvt. Ltd, Hyderabad, India. The raw sequence reads were screened using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for quality assurances and further analysed using FastP ver 0.20.1 (Chen et al., 2018) to remove the adapter sequences, low-quality reads (Q < 20) and reads with <30 bp length. After quality filtering, high-quality reads were assembled using animal mtDNA database sequences integrated in GetOrganelle ver 1.7.5 (Jin et al., 2020). The process involved following steps: (a) mapping reads to seed sequence against database using Bowtie 2 (Langmead & Salzberg, 2012), (b) conducting de novo assembly using SPAdes with different k-mer values (Bankevich et al., 2012), and (c) accurate identification of target organelle contigs. The assembled mitogenomes were validated by mapping against the red deer (Cervus elaphus) mitogenome (CerEla 1.0, OU343111.1).

Data analysis

Mitogenome organisation

A total of 215.84 million (female sample) and 284.22 million (male sample) sequence reads were generated through NGS approach. For mitogenome assembly, following paired and unpaired mapped reads were used: female-50,466 and 3,070; male-114,948 and 5,152. After data filtering, 213.7 million (female, 99.02%) and 281.78 million (male, 99.1%) of quality reads were obtained. The assembly lead to ~10× (female, 29.8 GB data) and ~13.74× (male, 39.3 GB data) coverage of Hangul mitogenome. The assembled mitogenome of Hangul generated in this study was of 16,351 bp length and matched with the available Hangul and other species reference sequences. The GenBank accession numbers for the sequences are ON416884 (female) and ON416885 (male), respectively.

The annotated Hangul mitogenome showed conserved gene order like other red deer species i.e. protein-coding genes (PCGs, n = 13), tRNA genes (n = 22), ribosomal genes (n = 2) and a non-coding control region (Table 2). Out of these, 29 genes are encoded on the heavy strand (H-strand) and remaining nine (NADH6 and eight tRNA genes: tRNA^Gln, tRNA^Asn, tRNA^Ala, tRNA^Cys, tRNA^Tyr, tRNA^Ser, tRNA^Glu and tRNA^Pro) are encoded on the light strand (L-strand). The number of overlapping regions (n = 13) and intergenic spaces (n = 15) were also consistent with other red deers (Table 2). NADH1, NADH2, NADH4L, NADH5, two cytochrome oxidases (COII and COIII), ATP8 and ATP6 and O_L, showed overlaps ranging from 1 to 40 bp. The tRNA^Leu gene showed the largest overlap (9 bp) with NADH5 whereas three tRNA genes (tRNA^Ile, tRNA^Thr, tRNA^Gly) showed the smallest overlap (1–3 bp) (Table 2). Similarly, the intergenic spaces ranged between 1–7 bp, with the longest between tRNA^Ser and tRNA^Asp, respectively. The O_L was 31 bp long and found in the WANCY region between tRNA^Asn and tRNA^Cys. Overall, the nucleotide composition comprised of 28.8% T, 24.5% C, 33.3% A and 13.5% G, with a positive AT value (0.072) (and a negative GC value of −0.289), indicating a AT-rich Hangul mitogenome (Table 1).

Protein coding genes

Hangul mitochondrial protein-coding genes (PCGs) was 11,354 bp in length with an AT-skewed base composition of 29.7% T, 25.9% C, 32.3% A and 12.2% G (Table 1). Out of the 13 coding genes majority (n = 12) were on the H-stand. These included six NADH genes (NADH1, NADH2, NADH3, NADH4, NADH5, NADH6 and NADH4L), ATPases (ATP6 and ATP8), cytochrome c oxidases (COI, COII and COIII) and cytochrome b (Cyt b). NADH4L was coded on the L-strand consistent with mitogenome organisation of other red deer species. NADH2, NADH3 and NADH5 have ATA start codon compared to others that have ATG. Only Cyt b (AGA), NADH2 (TAG), NADH4 and COIII had incomplete stop codon (T_ _) at the 5′ terminal of the adjacent gene to be completed during post transcription modification i.e. polyadenylation process. The Relative Synonymous Codon Usage (RSCU) for all PCGs consist of 3,784 codons.

rRNA, tRNA genes and control region

Hangul mitogenome has two rRNA genes: 12s rRNA (located between tRNA^Phe and tRNA^Val) and 16s rRNA (between tRNA^Val and tRNA^Leu) (Fig. 1). The total length of these genes is 2,526 bp with an average positive AT skew (0.218) and negative GC skew (−0.097) (Table 1). In addition, it also contains 22 tRNA genes (total length of 1,515 bp) with length varying from 60 bp (tRNA^Ser) to 75 bp (tRNA^Leu) (Table 1). Most of these tRNAs (n = 14) are located on H-strand and majority (n = 21) showed the typical secondary cloverleaf structures. Leucine and Serine have two tRNAs each and tRNA^Ser (GCT) lacks the D-paired arm. The Control region (CR) is 917 bp long and located between tRNA^Pro and tRNA^Phe. The base composition of CR was 31.7% T, 23.6% C, 29.4% A and 15.3% G with average negative AT (−0.038) and GC (−0.213) skew (Table 1).

Phylogenetic analysis and divergence time estimate

Mitogenome-based Bayesian phylogenetic reconstruction of the red deer complex (with 18 available species sequences including outgroups) reveals four mitochondrial lineages: C. elaphus, C. hanglu, C. canadensis and C. nippon. These lineages can be broadly divided into two major groups: the elaphoid clade (C. hanglu and C. elaphus) and the wapitoid clade (C. canadensis and C. nippon). The C. elaphus and C. hanglu forms the basal clade, followed by the C. canadensis and C. nippon as sister groups. The Hangul (C. hanglu hanglu) sequences (both generated in this study and the already available one) formed a distinct clade with its other subspecies C. hanglu yarkandensis and sister to the Hungarian red deer (C. elaphus hippelaphus). The tMRCA analysis suggested a divergence period spanning from 2.36 (Height Posterior Density (HPD) 4.23–1.06 million years ago (Mya)) to 0.16 (HPD 0.35–0.03 Mya) Mya (Fig. 2) for the red deer species complex. Our results indicated the divergence of the elaphoid group from the wapitoid group ~2.36 Mya (HPD 4.23–1.06 Mya), corresponding to the early Pleistocene period (Kuwayama & Ozawa, 2000). The Tarim deer species group separated ~1.55 Mya (HPD 2.88–0.58 Mya), overlapping with the Grünz glacial period in Europe (Kurten, 1968; Kuwayama & Ozawa, 2000). At Hangul group level, results suggest a coalescence of 0.75 Mya (HPD 1.51–0.18 Mya) between C. hanglu and C. yarkandensis.

The Hangul sequences generated from the male and female individuals were almost identical with five variable sites (n = 3 in PCGs and n = 2 in noncoding region). However, the Hangul mitogenome from GenBank showed 15 variable sites (n = 6 in the coding and n = 9 in the non-coding regions with three insertions in the rRNA region) with the sequences generated in this study. The mitogenomes of Hangul and closely related Cervus h. yarkandensis showed 345 variable sites. The highest difference was found with Cervus n. yakushimae (795 variable sites), but during annotation it was found that part of the mitogenome was missing and therefore was removed from further genetic distance analyses. While the protein coding regions among all the deer remained consistent (11,354 bp), major differences were found in the non-coding regions which varied from 913 bp for Cervus h. yarkandensis to 1,231 bp for Cervus n. centralis. One to seven base pair differences were seen for both tRNA (1,515–1,517 bp) and rRNA (2,525–2,532 bp) regions across the species (Table 1). As expected, the Hangul mitogenomes were genetically much closer (genetic distance = 0.001) than all other red deer species (genetic distance value varied from 0.022 to 0.052) (Table S1).