Sample collection and DNA extraction
We used 42 archival samples of northern red muntjac from a different region of India including northwest region (NWI=18), mainland (localities of north and central) (MLI=14), northeastern (NEI=3), and southern India (SI=5) and Andaman & Nicobar Islands (AI=2) (Table 3 and Fig. 8). Genomic DNA was extracted from tissue and hair samples using modified DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) protocol, whereas, for antlers and bone samples, we followed Gu-HCl based silica binding method [39]. These samples were collected from dead animals from the known localities by the local Forest Department of India and sent to Wildlife Forensic and Conservation Genetic Cell, Wildlife Institute of India, Dehradun. Since the samples were collected from the dead animals, Animal Ethics Committee approval was not required for this study. The authors confirmed that all experiments were performed following the relevant guidelines and regulations.
PCR amplification of complete mitogenome and sequencing
Polymerase chain reactions (PCRs) amplifications were carried out in 20μl volumes containing 10–20 ng of extracted genomic DNA containing, 1× PCR buffer (Applied Biosystem), 2.5mM MgCl2, 0.2 mM of each dNTP, 5 pmol of each primer, and 5 units of Taq DNA polymerase (Thermo Scientific). We performed PCR amplification using 23 overlapping fragments of complete mtDNA [40]. Besides, we included the fragments of complete Cytochrome b [41] and Cytochrome C oxidase-I gene [42] to increase the overlapping. PCR cycles for DNA amplification were 95°C for 5 min; followed by 35 cycles of 95°C for 40 sec. (denaturation), 54-56°C (annealing) for 40 sec, 72°C for 50 sec. (extension) and a final extension of 72°C for 15 min. The efficiency and reliability of PCR reactions were monitored by using control reactions. The PCR products were electrophoresed on 2% agarose gel and visualized under UV light in the presence of ethidium bromide dye. The amplified PCR products were treated with Exonuclease-I and Shrimp alkaine phosphatase (USB, Cleveland, OH) for 15 min. each at 37ºC and 80ºC, respectively, to remove any residual primer. The purified Amplicons were then sequenced bidirectionally using BigDye Terminator 3.1 on an automated Genetic Analyzer 3500XL (Applied Biosystems, Carlsbad, CA, USA).
Microsatellite loci amplification and genotyping
Nine cross-species microsatellite loci: INRA001 [43], Ca18 [44], BM6506 [45], RT1, RT27, NVHRT48 [46], CelJP27 [47], and T156, T193 [48] were selected for population genetic analysis of the red muntjac. Three sets of the multiplex set were carefully assembled based on molecular size and labeled fluorescent dyes of loci. Multiplex amplification was carried out in 10 μl reaction volumes containing 5 μl of QIAGEN Multiplex PCR Buffer Mix (QIAGEN Inc.), 0.2 μM labeled forward primer (Applied Biosystems), 0.2 μM unlabeled reverse primer, and 20–100 ng of the template DNA. PCR cycles for loci amplification were 95°C for 15 min; followed by 35 cycles of 95°C for 45 sec. (denaturation), 55°C (annealing) for 1 min, 72°C for 1 min. (extension) and a final extension of 60 °C for 30 min. Alleles were resolved in an ABI 3500XL Genetic Analyzer (Applied Biosystems) using the LIZ 500 Size Standard (Applied Biosystems) and analyzed using GeneMaker ver 2.7.4 [49].
Data Analysis
A total of 16 complete mitogenome sequences of red muntjac from five different localities of India were generated (Table 1). Sequences obtained from forward and reverse direction were edited and assembled using SEQUENCHER®version 4.9 (Gene Codes Corporation, Ann Arbor, MI, USA). The annotation of complete mitogenome was done using Mitos WebServer [50] and MitoFish [51]. Careful manual annotation was also conducted by sequence alignment with their related homologs sequences or species for ensuring the gene boundaries. Additionally, 36 sequences that consisted of northern (n=17); southern (n=17) and Sri Lankan (n=2) muntjac were included from GenBank to cover wide distribution ranges for better insight (Fig. 1; Supplementary table ST 1). The dataset of 52 sequences of red muntjac was aligned using the CLUSTAL X 1.8 multiple alignment programs [52] and alignments were checked by visual inspection. DnaSP v 5 [53] was used to estimate the haplotype (h) and nucleotide (π) diversity. A phylogenetic analysis was conducted in BEAST ver 1.7 [54]. For a better insight into the phylogenetic relationships, Bos javanicus (JN632606) was used as an out-group. The resulting phylogenetic trees were visualized in FigTree v1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/). The spatial distribution of haplotypes was visualized by a median-joining network, which was created using the PopART software [55]. The genetic distance between lineages was calculated based on the Tamura-3 parameter with a discrete Gamma distribution (TN92+G) with the lowest BIC score value implemented in MEGA X [56].
Estimating divergence dating
To estimate divergence times of red muntjac clades, we inferred genealogies using a relaxed lognormal clock model in BEAST ver 1.7 [57]. We performed the dating estimates using fifteen sequences downloaded from NCBI, including the species of Cervidae and Muntiacini. i.e., the Chital (A. axis, JN632599), Swamp deer (R. duvaucelii, NC020743), Red deer (C. elaphus, AB245427), European Roe deer (C. capreolus, KJ681491), Fallow deer (D. dama, JN632629), Water deer (H. inermis, NC011821), Mule deer (O. hemionus, JN632670), Hog deer (A. porcinus, MH443786), Formosan sambar (R.u.swinhoei, DQ989636 ), Tufted deer (E. cephalophus, DQ873526. ), Chinese muntjac (M. reevesi, NC_008491), Giant muntjac (M. vuquangensis, FJ705435), Putao muntjac (M. puhoatensis, MF737190) and Black muntjac (M. crinifrons, NC_004577) and the sequence of one species of the family Bovidae, i.e., the Banteng (B. javanicus, FJ997262). The TMRCA of bovids and cervids was set at a calibration point to 16.6 ± 2 million years ago (Mya) with a normal prior distribution [13, 58]. We used a Yule-type speciation model and the HKY + I + G substitution rate model for tree reconstruction. We conducted two independent analyses, using MCMC lengths of 10 million generations, logging every 1000 generations. All the runs were evaluated in Tracer v. 1.6. The final phylogenetic tree was visualized in FigTree v.1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).
Microsatellite analysis
A total of nine polymorphic loci were used to analyze the 42 samples of red muntjac for population genetic studies. The CERVUS ver 3.0.6 program [59] was used to estimate the polymorphic information content (PIC), the number of alleles per locus, the observed (Ho) and, expected (HE) heterozygosity. The allelic richness (Ar) and mean inbreeding coefficient (FIS) [60] was estimated using FSTAT ver 2.9.3 [61]. All the loci were checked for under HWE in GenAlEx v6.5 [62]. Factor correlation analysis (FCA) was performed using the GENETIX 4.05 software package [63]. CONVERT 1.31 [64] was used to convert the required input file format.
To determine the genetic structure of the population, the DAPC method was implemented using the ADEGENET package in R [65]. DAPC is a multivariate and model-free approach that maximizes the genetic differentiation between groups with unknown prior clusters, thus improving the discrimination of populations. This analysis does not require a population to be in HWE. The pairwise FST values (gene flow) among the populations were calculated using GenAlEx v6.5 [62].
Spatial genetic analysis
The correlation between the pairwise genetic and geographic distances was performed to detect the pattern of isolation by distance between the disjointed areas, according to Mantel’s test implemented in Alleles in Space 1.0 [66].