Skip to main content
Download PDF

Consideration of genetic variation and evolutionary history in future conservation of Indian one-horned rhinoceros (Rhinoceros unicornis)

BMC Ecology and Evolution volume 22, Article number: 92 (2022) Cite this article

Abstract

Background

The extant members of the Asian rhinos have experienced severe population and range declines since Pleistocene through a combination of natural and anthropogenic factors. The one-horned rhino is the only Asian species recovered from such conditions but most of the extant populations are reaching carrying capacity. India currently harbours ~ 83% of the global wild one-horned rhino populations distributed across seven protected areas. Recent assessments recommend reintroduction-based conservation approaches for the species, and implementation of such efforts would greatly benefit from detailed genetic assessments and evolutionary history of these populations. Using mitochondrial data, we investigated the phylogeography, divergence and demographic history of one-horned rhinos across its Indian range.

Results

We report the first complete mitogenome from all the extant Indian wild one-horned rhino populations (n = 16 individuals). Further, we identified all polymorphic sites and assessed rhino phylogeography (2531 bp mtDNA, n = 111 individuals) across India. Results showed 30 haplotypes distributed as three distinct genetic clades (Fst value 0.68–1) corresponding to the states of Assam (n = 28 haplotypes), West Bengal and Uttar Pradesh (both monomorphic). The reintroduced population of Uttar Pradesh showed maternal signatures of Chitwan National Park, Nepal. Mitochondrial phylogenomics suggests one-horned rhino diverged from its recent common ancestors ~ 950 Kya and different populations (Assam, West Bengal and Uttar Pradesh/Nepal) coalesce at ~ 190–50 Kya, corroborating with the paleobiogeography history of the Indian subcontinent. Further, the demography analyses indicated historical decline in female effective population size ~ 300–200 Kya followed by increasing trends during ~ 110–60 Kya.

Conclusion

The phylogeography and phylogenomic outcomes suggest recognition of three ‘Evolutionary Significant Units (ESUs)’ in Indian rhino. With ongoing genetic isolation of the current populations, future management efforts should focus on identifying genetically variable founder animals and consider periodic supplementation events while planning future rhino reintroduction programs in India. Such well-informed, multidisciplinary approach will be the only way to ensure evolutionary, ecological and demographic stability of the species across its range.

Peer Review reports

Background

The members of Rhinocerotidae family were once one of the most diverse and widely distributed terrestrial herbivores with complex evolutionary history [1]. By late Pleistocene, this family was reduced to only nine species (from more than 100 species) spread across Eurasia (seven species) and Africa (two species) [1, 2]. Subsequently, early Holocene global warming (after Last Glacial Maxima) triggered their extinction in western Eurasia and southward movement of eastern Eurasian rhinos, leading to their distribution across Southeast Asia [2, 3]. Further, the range of all Eurasian rhino species (Javan, Sumatran and One-horned rhino) were affected by a combination of natural and anthropogenic factors during Pleistocene-Holocene transition period [15–9 thousand years ago (Kya)] [1, 3,4,5,6], followed by recent events of exploitation of natural resources (during colonial era), industrialisation and poaching (since seventeenth century) [7,8,9,10]. Population size of the most widely distributed Javan rhinos (during Holocene) [11] were greatly reduced during human population expansion since 10,000 years ago [3], whereas the Sumatran rhino populations became fragmented and isolated (since Holocene) due to submerged Sundaland corridors (late Pleistocene) [6]. The one-horned rhinos faced climate-change driven habitat shrinkage in late Pleistocene [12]. Currently the Javan and Sumatran rhinos are categorized as Critically Endangered (~ 60 Javan rhino—[13] and < 100 Sumatran rhinos—[10]) and one-horned rhino as Vulnerable by IUCN (~ 3700 individual, [14]). Recovery of these species in their natural habitats requires deeper understanding of demography, ecology and genetics for appropriate conservation measures.

The one-horned rhino, being the only Asian species recovered from severe population decline in the past are critical for the evolutionary potential of this group. With a current population size of ~ 3700 individuals (increased from few hundred individuals in 1990s), it retains ~ 96% of the Asian rhino population [10, 13, 14]. As majority of the current one-horned rhino bearing areas in India and Nepal are reaching to their carrying capacities [15, 16], future conservation efforts are directing towards reintroduction-based programmes. Detailed genetic assessment of the existing rhino populations is critical in this regard since strong historical demographic declines has led to loss of genetic variation in all rhino species (Black rhino—[17], White rhino—[18], Sumatran rhino—[6], Javan rhino—[13]). For example, Liu et al. [1] suggested low population size and reduced genetic diversity across Rhinocerotidae family for an extended period of time. Similarly, mitogenome-based phylogeography reported low variation in both Sumatran [10] and Javan [13] rhinos, but no such data is available for one-horned rhinos.

In this paper, we investigated the phylogeography and evolutionary history of one-horned rhinos in India (henceforth Indian rhino) as it harbours 83% [19] of the global population of this species. We sequenced the polymorphic sites in the Indian rhino mitogenome in 111 wild individuals surveyed across seven extant populations covering the states of Assam, West Bengal and Uttar Pradesh. Further, we identified the Evolutionary Significant Units (ESUs) in Indian rhinos and suggested appropriate conservation measures to secure the evolutionary potential of this species. We believe that the results will provide the most exhaustive genetic information for Indian rhinos that would be useful in future reintroduction and population management efforts.

Results

Rhino mitogenome data and comparative analyses

Sequencing with 23 primers (Additional file 1: Table S1) generated 16,828 bp mitogenome (Additional file 2: Fig. S1) for wild Indian rhino (n = 16, Genbank: MZ736693–MZ736708, Additional file 1: Table S2). Comparison with the available one-horned rhino mitogenome data (Genbank: X97336) showed identical patterns of gene annotations. Composition analysis revealed AT-skewed mitogenome with 13 protein coding genes, 22 tRNA, 2 ribosomal genes and a non-coding control region (Additional file 1: Table S3). Comparative analyses with other rhino species (Additional file 1: Table S4 and Additional file 2: Fig. S2) revealed that the Indian rhinos have low segregating sites (SJava = 15,514, SAfrica = 10,680, SSumatra = 130, SIndia = 18) and nucleotide diversity (πJava = 0.56, πAfrica = 0.43, πSumatra = 0.003, πIndia = 0.0005) but high haplotype diversity (HdSumatra = 0.96, HdIndia = 0.93, HdJava = 0.91, HdAfrica = 0.67). Both African rhino species (white and black rhino) data were combined for this analyses as no intra-species variation was observed in the available data.

Phylogeography of wild Indian rhinos

Out of 15 primers designed to assess genetic variation, eight were finally used (Additional file 1: Table S1) to amplify all 21 polymorphic sites (covering 2531 bp sequence) of rhino mitogenome. This data was generated for additional 95 unique individuals (n = 56 tissue and 39 dung, Additional file 1: Table S2) (Genbank: MZ771364–MZ771458, MZ771459–MZ771553, MZ771554–MZ771648, MZ771649–MZ771743, MZ771744–MZ771838, MZ771839–MZ771933 and MZ771934–MZ772028). Remaining samples could not be used as they were rejected due to low amplification success for microsatellite data (n = 12), genetic recaptures (n = 13) and individuals from adjacent midden sites (n = 24). Sequencing comparison showed that out of the 21 polymorphic sites, two and three sites were specific to West Bengal and Uttar Pradesh, respectively, whereas all others were shared at different levels among the three states (shared among three states—10 sites, Assam–Uttar Pradesh—eight sites, Assam–West Bengal—four sites, West Bengal–Uttar Pradesh—0 sites, Additional file 1: Table S5). Median joining network (n = 111 individuals) showed a total of 30 haplotypes (h) across India. Majority of these haplotypes (93.3%, n = 28) were from Assam whereas both West Bengal (one haplotype, n = 20) and Uttar Pradesh (one haplotype, n = 10) populations were found to be monomorphic (Fig. 1). The sequence from Bihar rhino individual was identical to the Uttar Pradesh population. Population-wise genetic variation indices (Table 1) showed overall highest values for KNP (n = 46; S = 18, h = 19, π = 0.0021, Hd = 0.85), followed by MNP (n = 12; S = 14, h = 6, π = 0.0023, Hd = 0.89), ONP (n = 12; S = 9, h = 6, π = 0.0016, Hd = 0.89) and PWLS (10; S = 2, h = 3, π = 0.0002, Hd = 0.51). Bayesian genetic clustering corroborated with the earlier pattern (K = 3) where samples from West Bengal and Uttar Pradesh formed distinct clusters whereas Assam showed geographically intermixed fixed haplogroups (Fig. 1). The genetic differentiation (pairwise Fst) values among these three clusters were significantly high ranging from 0.68 to 1 (Table 2, indicating highly structured populations). The hierarchical AMOVA analysis using two separate groupings: (a) seven populations and (b) three states showed higher within population (50%) and between group variance (45%) (Table 2). Such pattern indicates that overall genetic structure is influenced by differentiation at clade level and the amount of diversity present within the Assam clade (Tables 1 and 2).

Fig. 1
figure 1

Representation of mtDNA variations and genetic structure in Indian rhinos based on 2531 bp concatenated sequence covering all polymorphic sites across seven genes. a Median joining network with park-level colour codes; b Haplotype frequencies at each of the sampled areas covering all variations (n = 30 haplotypes); c Bayesian clustering shows monomorphism in Uttar Pradesh (with sample from Bihar, n = 11) and West Bengal (n = 20) populations and polymorphism in Assam (n = 80)

Full size image
Table 1 mtDNA diversity indices of all rhino populations in India (n = 111)
Full size table
Table 2 Results of pairwise genetic differentiation and hierarchical AMOVA test (Bihar sample considered under Uttar Pradesh clade)
Full size table

Divergence time of different rhino clades and demographic history

The Bayesian phylogeny showed similar pattern of three clades consisting of West Bengal, Assam (nodes C–E) and Uttar Pradesh (along with the Bihar sample, Fig. 2). Based on the calibrated root nodes and Indian rhino-specific mutation rate (1.2 × 10–4 mean rate of substitution per site per million years, Additional file 2: Fig. S3), tMRCA analysis suggested a divergence period spanning from 950 (HPD 1360–810 Kya) to 50 (150–10 Kya) Kya (Fig. 2). Our results indicated the divergence of Indian rhinos ~ 950 Kya (node A, Fig. 2) corresponds to the emergence period of one-horned rhino ancestors in the subcontinent [5, 12]. Next, the Assam population diverged from the remaining clades at ~ 500 Kya (HPD 680–330 Kya, nodes B & C, Fig. 2). This is supported by reports of multiple rhino movements away from Assam (along Siwalik as well to Siva-Malayan region) during this period [12, 20]. At population level, results suggest a relatively earlier coalescence of Assam ~ 190 Kya (HPD 300–70 Kya, node D & E, Fig. 2) compared to West Bengal and Uttar Pradesh (~ 50 Kya, HPD 150–10 Mya, node F & G, Fig. 2). This period (120–10 Kya) is known for confinement of rhinoceros to the north and north-east of India due to monsoon intensification and grassland dominance [5, 12, 21].

Fig. 2
figure 2

Phylogenetic relationship and assessment of divergence time in Indian rhino populations. The left pane shows the clustering of three maternal clades of West Bengal samples (green), Uttar Pradesh (blue) and Assam (red). Javan rhino sequence was used as outgroup. The posterior probability values (≥ 0.9) are shown in bold. The right pane indicates the divergence of Indian rhinos ~ 0.95 Mya, where the Assam population coalesce first (~ 0.19 Mya), followed by divergence of West Bengal and Uttar Pradesh (0.06–0.05 Mya). Node-specific ages are marked (with posterior probability values ≥ 0.9). The major corroborating paleobiogeographical events are presented above

Full size image

All four BSP analyses showed similar population trends for overlapping time periods where the combined data identified trends at a deeper coalescence period compared to the clade-specific data (Fig. 3). The Assam clade showed a steep increase in female effective population size ~ 110 Kya followed by constant population size from ~ 90 Kya (Fig. 3a) whereas West Bengal and Uttar Pradesh clades showed similar demographic trends of stable populations from ~ 60 Kya (Fig. 3b, c). The combined dataset showed a steep decline in population size ~ 300–200 Kya followed by a gradual increase ~ 120 Kya and steep rise ~ 60 Kya (Fig. 3d).

Fig. 3
figure 3

Bayesian skyline plot analysis (BSP) to determine the changes in female effective population size across three clades, a Assam, b West Bengal, c Uttar Pradesh and d combined dataset of Indian rhinos. The vertical lines represent the HPD intervals of the given divergence time for each analysis whereas the shaded horizontal area is the HPD of the median effective size value

Full size image

Discussion

This study presents the most extensive mitochondrial DNA phylogeography of one-horned rhinos across its Indian distribution. Careful considerations involving mitogenome sequencing of representative samples across Indian rhino-bearing areas, identification of all polymorphic regions and their amplification from spatially-covered rhino samples helped us achieving accurate assessment of mtDNA variations. To the best of our knowledge, this is the first report of wild Indian one-horned rhino mitogenome from all the extant populations. Despite relatively similar haplotype diversity of Asian rhinos [India—0.93 (16 samples), Sumatra—0.96 (15 samples), and Java—0.91 (6 samples), respectively], Indian rhino mitogenome showed much lower values for segregating sites and nucleotide diversity (Additional file 1: Table S4). Such mitogenome comparisons may be affected by limited sample size (earlier studies in African rhinos have reported higher diversity based on partial mitogenome data with more samples [17, 18]) or representation of historical genetic variations (in Javan rhinos, [13]). However, it was surprising to observe that despite similar historical demographic incidences (severe population decline due to habitat shrinkage [6, 12] and anthropogenic pressures [9, 10]) Indian rhino retain much lower genetic variation than their Sumatran counterpart. This can be potentially attributed to recovery of the Indian species from extremely low founder population (as indicated by high Hd but low π) [22, 23].

As expected, the phylogeography data (2531 bp mtDNA, n = 16 samples) revealed higher number of haplotypes than the mitogenome data (n = 30 haplotypes) (due to large sample size). The only other study on one-horned rhino mtDNA variations (based on partial control region sequences, 428 bp) reported 10 haplotypes (Kaziranga National Park, India—4 and Chitwan National Park, Nepal—6, respectively) and moderate level of genetic difference (Fst value of 0.39 between them) [24]. Careful scrutiny of our data revealed that all the polymorphic sites (or identified segregating sites) were found in fixed positions within one-horned rhino mitogenome (Additional file 1: Table S5) across India. Given the distribution of polymorphic sites in the sequenced mitogenomes and our sampling coverage, it is likely that these data represent the majority or perhaps all extant mtDNA haplotypes in Indian rhinoceros populations. This claim is also supported by the similar haplotype diversity values from the mitogenome and the phylogeography datasets (0.93 and 0.9, respectively). Our study also shows that the Indian rhinos have the highest number of haplotypes compared to the other genus/species reported so far [10, 13, 17, 18]. The clustering analysis of the concatenated rhino sequences showed three distinct genetic clades (corresponding to the states of Assam, West Bengal and Uttar Pradesh) with high Fst value (0.68–1), corroborating with the haplotype network patterns. Mantel test (− 0.83, p = 1) confirmed that such strong genetic structuring is not due to isolation by distance pattern, but driven by lineage-specific evolutionary history (as suggested by AMOVA results). Such pattern of higher within population and between group variance (50% and 45% in Indian rhinos, respectively) indicates that the mitochondrial genetic variation observed in extant Indian rhino is influenced by both evolutionary diversification and retention of diversity at population level only for Assam clade. As two of the clades are monomorphic, they contribute very less proportion of among-population within-group variations (5%). Similar data has also been described in other species such as barking deer—[25], dog—[26] etc. Interestingly, we found that the sequence from the Bihar sample (representing samples from Nepal) was identical to the Uttar Pradesh sequences, including the state-specific SNPs. This pattern was expected as the founder animals of the reintroduced Uttar Pradesh population were sourced from Chitwan National Park, Nepal (four dominant breeding females) and Pobitora Wildlife Sanctuary of Assam (dominant breeding male) [27]. Further comparison of 13 partial D-loop sequences from Chitwan National Park, Nepal Zschokke et al. [24] confirmed this pattern, indicating that the mtDNA signature of the Uttar Pradesh population belongs to Nepal. Given that the entire Uttar Pradesh rhino population showed only one haplotype, future studies need to evaluate the mtDNA variation in the Nepal population.

The phylogenetic analyses reconfirmed the relationship among the existing members of the Rhinocerotidae family [10, 28, 29] where the Sumatran and African rhino formed sister clades, separated from the Rhinoceros sp. based on the extant rhino genus/species sequence data only (Woolly rhinoceros sequence was not used). The within species tree topology corroborated with the haplotype network results as Assam and Uttar Pradesh formed phylogenetically closer clades as compared to West Bengal. We believe that the observed phylogenetic pattern of West Bengal being separate clade is influenced by lesser shared polymorphic sites between West Bengal and other two clades (Additional file 1: Table S5). Combined together, we interpret that the one-horned rhino diverged from its recent common ancestors ~ 950 Kya and different populations (Assam, West Bengal and Uttar Pradesh/Nepal) coalesce around ~ 190–50 Kya time period (Fig. 2). The molecular dates were comparable to other published literature on rhino evolution [5, 12, 21] and supported by the paleobiogeographic history of the Indian subcontinent [12, 21]. For instance, the inward movement of rhinos from Assam along Siwalik (680–330 Kya, node B & C) coincides with drop in the sea level which facilitated movement of multiple genera (for example, Elaphas, Panthera, Rhinoceros, Muntiacus etc.) through Siva-Malayan route [20, 30]. Report of one-horned and Javan rhino co-existence in Bhutan ~ 560 Kya [13] provide further support of such movements. Finally, the coalescence time of the three Indian clades (~ 190–50 Kya, Fig. 2) corresponds to Holocene climatic optimum period known for monsoon intensification in north and north-east part of India resulting in range contraction for grassland dependent species [5, 12, 31]. We feel that our approach of using taxon-specific mutation rate and fossil data for node calibration has resulted in achieving such meaningful estimates of tMRCA. Future efforts should try to include molecular data from historical/ancient samples to tighten the variance associated with divergence estimates [32]. Overall, this approach reiterates the critical importance of large datasets (whole mitogenome from multiple individuals in this case), informative prior settings and its assessment with posterior outputs, taxon-specific mutation rate, node calibration points etc. for accurate tMRCA estimation [33,34,35,36,37].

The BSP results with different datasets (combined vs. three individual clades) showed similar patterns of changes at different evolutionary timescale. The combined data indicate a historical decline in maternal effective population size ~ 300–200 Kya, followed by increasing tends during ~ 110–60 Kya (coinciding with Holocene climatic optimum period, also seen in the Assam clade analysis, Fig. 3). This pattern is similar to earlier findings (described in [1] based on whole genome data) with a difference in the Ne values arising from lower effective size in mtDNA [38]. The West Bengal and the Uttar Pradesh clades did not show any changes in population trajectories owing to the monomorphic data. It is noteworthy to point out that such mitochondrial DNA-based analyses would only capture the demographic history at longer evolutionary time scale, and use of suitable nuclear markers (microsatellites, SNPs etc.) could provide much powerful demographic inferences [38].

The spatially exhaustive sampling coverage and the patterns of population structure brings out some critical conservation perspectives for the Indian rhinos. The phylogeographic and mitophylogenomic patterns suggest three distinct clades with state-specific evolutionary histories. As these populations are morphologically undistinguishable and interbreed among themselves (Dudhwa, Uttar Pradesh population is genetically mixed, [27]), we suggest that they should be recognised as ‘Evolutionary Significant Units (ESUs)’ [39]. It is therefore important to use such information towards conservation and management efforts for each of these populations [39,40,41]. Despite strong recoveries across all existing populations since late 1990s, recent analyses suggest high extinction probability of the species [42], and further conservation efforts are mostly concentrated on translocation activities [16, 43]. Till date genetic information of the species has not been used in translocation planning (possibly due to lack of sufficient data), and the genetic signatures described in this study would be very helpful to increase variation in target populations. For example, the Uttar Pradesh and West Bengal population show state-specific monomorphic haplotypes representing unique but genetically depauperate populations. Based on the data presented here, suitable founder animals from Assam populations can be considered for future translocation programs in these areas, thereby increase the genetic diversity of these populations to combat any potential stochastic events [40, 41]. However, such efforts would impact the suggested ESU categorizations due to mixing of different gene pools among populations. Another important aspect for management consideration would be better planning for translocation events to any of the existing or new areas [16, 43]. For example, the reintroduced rhino population in Assam (Manas National Park) showed much higher mtDNA variation (six haplotypes), possibly due to periodic supplementation of individuals of varied genetic ancestry across different wild rhino populations [44] compared to Dudhwa NP (single haplotype) of Uttar Pradesh (single supplementation event). As multiple reintroduction programs are planned as per the ‘National Conservation Strategy for the Indian One horned rhinoceros (Rhinoceros unicornis), Government of India, Ministry of Environment Forest and Climate Change, 2021’ objectives (in the states of Uttar Pradesh, Bihar, West Bengal and Assam) in near future, we suggest that all future efforts should adopt the Manas NP model with consideration of selecting genetically variable founder animals, multiple reintroduction events etc. [16, 45].

Conclusion

The one-horned rhino was found throughout the Indo-Gangetic plains during the early twentieth century [43, 46] but faced drastic reductions in distribution and population size (including local extinctions) [47, 48], followed by one of the most successful species recovery (increase in population size) in wild across the world [7, 43, 48]. We present the first assessment of range-wide mitogenome diversity in Indian rhinos where we emphasize the importance of large data, spatial sampling coverage of populations and evolutionary history as fundamental information for future population reintroduction/recovery programs (as suggested in case of other species [49,50,51,52]). Our results are important for Indian rhino conservation because they suggest higher genetic diversity than earlier reported [24]. However, the existing habitats are small, disjunct, isolated and reaching their respective carrying capacities [16, 48] and conservation options are becoming limited except establishing new habitats and translocation-driven population enhancement [16]. We believe that the genetic information provided here will assist in identifying appropriate source populations and maintain adequate genetic diversity in the existing (and new) rhino populations, thereby ensuring evolutionary, ecological and demographic stability for their future survival.

Methods

Permission and ethical considerations

Data generated in this study is part of a collaborative programme titled “Implementing Rhino DNA Indexing System to counter rhino poaching threat and aid population management in India” (henceforth RhoDIS-India). Biological sampling from all the three rhino bearing states was permitted by Ministry of Environment, Forests and Climate Change (MoEF&CC), Government of India (Letter No. 4-22/2015/WL). Permission for dung sampling was provided by state forest departments of Assam (Letter No. A/GWL/RhoDIS/2017/913, 3653/WL/2W-525/2018, WL/FE.15/22), West Bengal (Letter No. 3967/WI/2W-525/2018) and Uttar Pradesh [Letter No. 1978/23-2-12 (G)]. We have also received one tissue sample from Valmiki National Park (henceforth NP), Bihar forest department assumed to be representing the wild rhinos of Nepal (Letter/no.-1296 dated 16.10.2020). No ethical permissions were required for tissues as they were collected from naturally dead rhinos as well as for dung samples.

Study area

During the 1600s, the one-horned rhinos were distributed throughout the northern Indian subcontinent covering all the major river basins from Pakistan to Indo-Myanmar borders. The species has lost most of its habitat and population size due to a range of anthropogenic interventions (habitat loss, hunting, poaching etc.) [7, 9, 14] and are currently distributed across only 12 protected areas covering > 2000 km2 area in India and Nepal [14]. This study on the Indian rhino was conducted across all extant rhino-bearing parks (n = 7) found across the states of Assam (n = 4 parks), West Bengal (n = 2 parks) and Uttar Pradesh (n = 1 park) situated in Terai-Duars region of north and north-east India (Fig. 4). Assam currently hosts the largest population of Indian rhinos (~ 80% of the total population), which are found across four isolated parks: Kaziranga NP, Orang NP, Pobitora Wildlife Sanctuary (WLS) and Manas NP. All of these populations have experienced severe hunting and poaching pressures during early-late 1990s, leading to local population sizes ranging from 50-few hundreds and local extinction in Manas NP, but have revived to their current population sizes [7, 14]. Between 2010 and 2020 Manas received 35 translocated rhinos from Kaziranga NP and Pobitora WLS part of the population recovery program [44]. The most recent population estimates of these parks are as follows: Kaziranga NP: ~ 2500; Orang NP: 100; Pobitora WLS: 101; and Manas NP: 42 [14].

Fig. 4
figure 4

Map of the study area and distribution of the final samples used in this study (n = 111). The top plate shows the position of the rhino-bearing parks across three Indian states (Uttar Pradesh, West Bengal and Assam). The reference sample of wild rhino received from the state of Bihar (*) is also presented

Full size image

The state of West Bengal currently retains ~ 350 rhino individuals distributed between two parks, Gorumara NP (52 individuals) and Jaldapara NP (> 250 individuals). This population has recovered from a severe population decline of ~ 20 individuals during early 1900s (due to severe habitat loss) [9]. The rhino population in Uttar Pradesh was locally extinct along with the entire Terai in mid 1990s (mostly due to habitat loss and hunting). During 1984–85, rhinos were reintroduced in Dudhwa NP (Uttar Pradesh) from Chitwan NP, Nepal and Pobitora WLS, Assam and currently this park hosts ~ 40 rhinos [27]. Apart from this, wild rhinos are occasionally reported from Valmiki NP, Bihar (adjacent to Chitwan NP, Nepal in the Indian part of Terai). These rhinos are either swept down by natural flooding or use the grasslands along the river Gandak within Valmiki NP during monsoon seasons.

Biological sampling

Overall the sampling strategy in this phylogeography study was to select unique rhino individuals from different parts of the species distribution in India. A total of 160 samples (72 tissues and 88 dung) covering four states (Fig. 4) were used to assess rhino mitochondrial genetic diversity. The tissue samples of naturally dead rhino were provided by respective forest departments as part of RhoDIS-India protocol (2017–2021). Further dung collection was done to ensure spatial coverage for areas with no representative tissue samples. Rhino dung sampling can be challenging in the wild due to their use of communal latrine system (middens) [53, 54]. In this study, sampling was conducted by intensive foot and vehicle surveys from already known midden sites across six rhino bearing parks (except Kaziranga NP). During sampling, only the fresh bolus from top of the midden was selected and swabbed twice with separate PBS-soaked sterile cotton swabs (Himedia, Mumbai, India). All samples were geo-tagged and transferred to laboratory in − 20 °C freezer till downstream processing.

DNA extraction

Tissue DNA was extracted using already established protocol for Indian rhino mentioned in Ghosh et al. [55]. For dung samples, a modified protocol from Biswas et al. [56] was used. In brief, samples were digested overnight with a combination of 700 μl ATL and 65 μl Proteinase K (20 mg/ml) at 56 °C, followed by QIAamp DNA Tissue Kit (QIAGEN Inc., Hilden, Germany) protocol with adjusted volumes. DNA was eluted twice in 100 μl preheated (70 °C) 1X TE buffer and stored in − 20 °C freezer. Extraction negative was used for each set of extraction (n = 23) to monitor possible contamination.

PCR amplification and sequencing

To assess genetic variation of the extant rhino populations, complete mitogenome data was generated for representative samples from each park (n = 15, see Additional file 1: Table S2 for details) and one from the Valmiki National Park, Bihar. These samples were selected based on their geographic locations representing the farthest samples within each park to ensure inclusion of potentially unrelated individuals. Mitogenome sequencing was performed using already published 23 overlapping primers [57]. For annealing temperature standardisation, gradient PCR was set in 10 μl reactions containing 4 μl of 2× Qiagen PCR buffer mix (QIAGEN Inc., Hilden, Germany), 1 μl of primer (3 μM), 2 μM BSA (4 mg/ml), 1.4 μl of RNase free water and 5 ng of rhino tissue DNA. PCR conditions included an initial denaturation (95 °C for 15 min); 35 cycles of denaturation (95 °C for 30 s), annealing (50–60 °C gradients for 40 s) and extension (72 °C for 40 s); followed by a final extension (72 °C for 10 min). During each set of reactions, PCR and extraction negatives were included to monitor contamination. Amplified products were visualized with 2% agarose gel, cleaned with Exonuclease (Thermo Scientific, Waltham, USA) and Shrimp Alkaline Phosphatase (Amresco, Solon, USA) mixture and sequenced bidirectionally in an ABI 3500XL bioanalyzer (Applied Biosystems). Out of these 23 primers, two did not show amplification in any samples. The remaining sequences (n = 21 from 16 individuals) were aligned with the available one-horned rhino mitogenome (Genbank: X97336, [28]) in Mega v7 [58]. Two primers were designed manually in the flanking conserved regions adjacent to the gaps (Additional file 1: Table S1) and sequences were generated from all the samples (n = 16).

The complete mitochondrial sequences (n = 16) were aligned and manually screened to identify the segregating sites. Further, a total of 15 primers were designed (multiple primers covering the segregating sites) to amplify all the polymorphic sites as < 500 bp fragments to ensure higher success rate from—poor quality dung DNA samples. These primers were standardised following same protocol described above. For all field collected samples (tissue = 56 and dung = 88) individual identification was performed using a panel of 14 microsatellites (described in [55]). After PCR amplification and genotyping of the markers, samples with 12–14 loci data were selected for downstream analysis and genetic recaptures were removed. To ensure removal of closely related individuals in our dataset we selected one sample from adjacent midden sites. Sequence data (2531 bp covering seven genes) was generated for the selected individuals to assess phylogeography patterns.

Complete mitogenome annotation and comparative analysis

All rhino sequences (n = 16) were aligned in Mega v7 to generate a complete mitogenome sequence and manually checked to identify any nucleotide ambiguities. Annotation was done using MITOS2 web with default settings and vertebrate mitochondrial genetic code [59] followed by mitogenome map construction with OGDRAW [60]. The mitogenome annotation was further confirmed with earlier published one-horned rhino mitogenome data (Genbank: X97336, [28]). To ascertain species-wise mitochondrial DNA diversity these sequences were aligned with already published rhino mitogenome sequences from Diceros bicornis (n = 2, Genbank: FJ905814, NC012682 [61]), Ceratotherium simum (n = 2, Genbank: Y07726, NC001808 [62]), Dicerorhinus sumatrensis (n = 15, Genbank: MF066629-MFO66643 [10] and Rhinoceros sondaicus (n = 6, Genbank: FJ905815 [61], MK909142, MK909146, MK909148, MK909149, MK909151 [13]). We calculated number of segregating sites (S), nucleotide (π) and haplotype diversity (Hd) using DnaSP v.5 [63] for all genes in the mitogenome.

Genetic diversity in Indian rhinos

Population-wise basic indices of genetic variations (S, π and Hd) were calculated for concatenated sequence data (2531 bp from seven genes) using DnaSP v.5 followed by a median joining [64] haplotype network constructed in PopART v. 1.7 [65]. To ascertain any possible population structure a Bayesian approach implemented in BAPS v.5.3 was used as it considers linked loci data [66]. Pairwise Fst and differential hierarchical AMOVA analysis was performed using Arlequin v. 3.0 [67] to confirm the pattern found in BAPS analysis.

Estimation of clade-specific divergence times and demographic history

To identify the clades, Bayesian phylogeny was constructed with MrBayes v. 3.2.7 [68] using 16 Indian rhino mitogenome and Javan rhino sequence (outgroup, as they are the sister clade of one-horned rhinoceros) [69]. Analysis was conducted using GTR + G substitution model determined by jModelTest v2.1.3 [70] (based on Akaike Information Criteria). The MCMC parameters included 2 runs of four chains each of 15 million generations with sampling after 1000 generations till split frequencies were below 0.01. Posterior probabilities were calculated for each node.

To estimate divergence among clades, rate of mutation for Indian rhino was calculated using BEAST v.2.3.6 [32]. Analysis was performed with five extant rhino mitogenome (without D-loop) (n = 11 sequences, India = 7 (haplotypes representing maximum variation in the data), Java = 1, Sumatra = 1, White = 1, Black = 1) along with horse (Equus caballus, Genbank: NC001640), donkey (Equus asinus, Genbank: NC001788), Asiatic wild ass (Equus hemionus, Genbank: NC016061) and zebra (Equus zebra, Genbank: NC018780) as outgroups. GTR + G substitution model was selected through jModelTest v2.1.3 for this multi-species data. Birth–death speciation was considered as tree prior [10, 34] along with uncorrelated relaxed log normal clock [10, 33]. During analysis, four established internal node calibration points (based on fossil records) with normal distribution priors were employed: (i) Caballine split (4 ± 0.5 million years ago (Mya)) [71, 72]; (ii) late Oligocene diversification of rhino groups (26 ± 3.5 Mya) [73]; (iii) split of rhinoceros genus (3 ± 0.5 Mya) [1]; (iv) origin of the perissodactyls (55 ± 3 Mya) [74, 75]. The first three calibration points were considered as monophyletic constraint [33] as the last point includes both ingroup and outgroup taxa.

tMRCA (time to Most Recent Common Ancestor) was inferred using the estimated mutation rate with lognormal distribution under strict molecular clock (intra species data, n = 16) [34, 35]. MCMC runs included 100 million generations, sampled at every 10,000 states with 10% burn-in. Data convergence was checked with Tracer v. 1.5 [76] and the final tree (with maximum clade credibility) was estimated with TreeAnnotator [77] and visualised using FigTree v.1.4.2 [78].

To estimate past fluctuations in population size, Bayesian skyline analysis (Bayesian skyline plot or BSP) was conducted using concatenated sequence data with monophyletic constraint to the identified maternal clades ensuring phylogenetic construction. Analysis was conducted with multiple datasets (each clade and combined data, respectively) to ascertain any possible impacts of genetic structure in the data [38]. In all cases coalescent BSP tree prior was used along with strict molecular clock, estimated mutation rate and clade specific divergence date. MCMC parameter settings and data convergence were identical to the tMRCA analysis.

Availability of data and materials

The dataset generated in this study is available in GenBank with Accession numbers MZ736693–MZ736708 (whole mitogenome sequences) and MZ771364–MZ771458, MZ771459–MZ771553, MZ771554–MZ771648, MZ771649–MZ771743, MZ771744–MZ771838, MZ771839–MZ771933 and MZ771934–MZ772028 (fragmented sequence data for phylogeography). The additional figures and tables are provided as Additional file 1: Tables S1–S4 and Additional file 2: Figs. S1–S3).

Abbreviations

BSP:

Bayesian skyline plot

ESUs:

Evolutionary significant units

Kya:

Thousand years ago

mtDNA:

Mitochondrial DNA

Mya:

Million years ago

Ne:

Effective population size

NP:

National Park

tMRCA:

Time to the most common recent ancestor

References

  1. Liu S, Westbury MV, Dussex N, Mitchell KJ, Sinding MHS, Heintzman PD, Duchêne DA, Kapp JD, von Seth J, Heiniger H, Sánchez-Barreiro F. Ancient and modern genomes unravel the evolutionary history of the rhinoceros family. Cell. 2021;184:1–12. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cell.2021.07.032.

    Article  CAS  Google Scholar 

  2. Wan X, Zhang Z. Climate warming and humans played different roles in triggering Late Quaternary extinctions in east and west Eurasia. Proc R Soc B. 2017;284:20162438. https://0-doi-org.brum.beds.ac.uk/10.1098/rspb.2016.2438.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Li X, Jiang G, Tian H, Xu L, Yan C, Wang Z, Wei F, Zhang Z. Human impact and climate cooling caused range contraction of large mammals in China over the past two millennia. Ecography. 2015;38:74–82. https://0-doi-org.brum.beds.ac.uk/10.1111/ecog.00795.

    Article  Google Scholar 

  4. Alroy J. A multispecies overkill simulation of the end-Pleistocene megafaunal mass extinction. Science. 2001;292:1893–6. https://0-doi-org.brum.beds.ac.uk/10.1126/science.1059342.

    Article  CAS  PubMed  Google Scholar 

  5. Lal P. Indica a deep natural history of the Indian subcontinent. In: The carnival of mammals. Random House India Pvt. Ltd: Delhi; 2016. p. 284–315.

    Google Scholar 

  6. Mays HL Jr, Hung CM, Shaner PJ, Denvir J, Justice M, Yang SF, Roth TL, Oehler DA, Fan J, Rekulapally S, Primerano DA. Genomic analysis of demographic history and ecological niche modeling in the endangered Sumatran rhinoceros Dicerorhinus sumatrensis. Curr Biol. 2018;28:70–6. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cub.2017.11.021.

    Article  CAS  PubMed  Google Scholar 

  7. Menon V. Under Siege: poaching and protection of greater one-horned rhinoceroses in India, TRAFFIC India; 1996. https://portals.iucn.org/library/fr/node/7008.

  8. Fernando P, Polet G, Foead N, Ng LS, Pastorini J, Melnick DJ. Genetic diversity, phylogeny and conservation of the Javan rhinoceros (Rhinoceros sondaicus). Conserv Genet. 2006;7:439–48. https://0-doi-org.brum.beds.ac.uk/10.1007/s10592-006-9139-4.

    Article  CAS  Google Scholar 

  9. Das PK, Borthakur U, Sarma HK, Talukdar BK. Population genetic assessment of extant populations of greater one-horned rhinoceros (Rhinoceros unicornis) in India. Eur J Wildl Res. 2015;61:841–51. https://0-doi-org.brum.beds.ac.uk/10.1007/s10344-015-0960-2.

    Article  Google Scholar 

  10. Steiner CC, Houck ML, Ryder OA. Genetic variation of complete mitochondrial genome sequences of the Sumatran rhinoceros (Dicerorhinus sumatrensis). Conserv Genet. 2017;19:397–408. https://0-doi-org.brum.beds.ac.uk/10.1007/s10592-017-1011-1.

    Article  CAS  Google Scholar 

  11. Antoine PO. Pleistocene and Holocene rhinocerotids (Mammalia, Perissodactyla) from the Indochinese Peninsula. C R Palevol. 2012;11:159–68. https://0-doi-org.brum.beds.ac.uk/10.1016/j.crpv.2011.03.002.

    Article  Google Scholar 

  12. Patnaik R. Neogene-Quaternary Mammalian paleobiogeography of the Indian subcontinent: an appraisal. C R Palevol. 2016;15:889–902. https://0-doi-org.brum.beds.ac.uk/10.1016/j.crpv.2015.11.004.

    Article  Google Scholar 

  13. Margaryan A, Sinding MHS, Liu S, Vieira FG, Chan YL, Nathan SKSS, Moodley Y, Bruford MW, Gilbert MTP. Recent mitochondrial lineage extinction in the critically endangered Javan rhinoceros. Zool J Linn Soc. 2020;190:372–83. https://0-doi-org.brum.beds.ac.uk/10.1093/zoolinnean/zlaa004.

    Article  Google Scholar 

  14. Talukdar B. Asian rhino specialist group chair report/Rapport du Groupe de Spécialistes du Rhinocéros d’Asie. Pachyderm. 2021;62:29–34.

    Google Scholar 

  15. Subedi N, Lamichhane BR, Amin R, Jnawali SR, Jhala YV. Demography and viability of the largest population of greater one-horned rhinoceros in Nepal. Glob Ecol Conserv. 2017;12:241–52. https://0-doi-org.brum.beds.ac.uk/10.1016/j.gecco.2017.11.008.

    Article  Google Scholar 

  16. Jhala HY, Qureshi Q, Jhala YV, Black SA. Feasibility of reintroducing grassland megaherbivores, the greater one-horned rhinoceros, and swamp buffalo within their historic global range. Sci Rep. 2021;11:1–15. https://0-doi-org.brum.beds.ac.uk/10.1038/s41598-021-83174-4.

    Article  CAS  Google Scholar 

  17. Moodley Y, Russo IRM, Dalton DL, Kotzé A, Muya S, Haubensak P, Bálint B, Munimanda GK, Deimel C, Setzer A, Dicks K. Extinctions, genetic erosion and conservation options for the black rhinoceros (Diceros bicornis). Sci Rep. 2017;7:41417. https://0-doi-org.brum.beds.ac.uk/10.1038/srep41417.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Moodley Y, Russo IRM, Dalton DL, Kotzé A, Smith S, Stejskal J, Ryder OA, Hermes R, Walzer C, Bruford MW. Contrasting evolutionary history, anthropogenic declines and genetic contact in the northern and southern white rhinoceros (Ceratotherium simum). Proc R Soc B. 2018;285:20181567. https://0-doi-org.brum.beds.ac.uk/10.1098/rspb.2018.1567.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Kakati P, Paine SK, Bhattacharjee CK, Bhattacharyya C, Sharma A, Phukan D, Basu A. Gut microbiome architecture of wild greater one-horned rhinoceros: a vulnerable species from Kaziranga National Park. India J Genet. 2021;100(2):1–7. https://0-doi-org.brum.beds.ac.uk/10.1007/s12041-021-01326-x.

    Article  CAS  Google Scholar 

  20. van den Bergh GD, de Vos J, Sondaar PY. The Late Quaternary palaeogeography of mammal evolution in the Indonesian Archipelago. Palaeogeogr Palaeoclimatol Palaeoecol. 2001;171:385–408. https://0-doi-org.brum.beds.ac.uk/10.1016/S0031-0182(01)00255-3.

    Article  Google Scholar 

  21. Patnaik R. Indian neogene Siwalik mammalian biostratigraphy. In: Fossil mammals of Asia. New York: Columbia University Press; 2013. p. 423–44.

    Chapter  Google Scholar 

  22. Avise JC. Phylogeography: the history and formation of species. London: Harvard University Press; 2000.

    Book  Google Scholar 

  23. da Silva M, Noll FB, et al. Phylogeographic analysis reveals high genetic structure with uniform phenotypes in the paper wasp Protonectarina sylveirae (Hymenoptera: Vespidae). PLoS ONE. 2018;13: e0194424. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0194424.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Zschokke S, Armbruster GF, Ursenbacher S, Baur B. Genetic differences between the two remaining wild populations of the endangered Indian rhinoceros (Rhinoceros unicornis). Biol Conserv. 2011;144:2702–9. https://0-doi-org.brum.beds.ac.uk/10.1016/j.biocon.2011.07.031.

    Article  Google Scholar 

  25. Martins RF, Fickel J, Le M, Van Nguyen T, Nguyen HM, Timmins R, Gan HM, Rovie-Ryan JJ, Lenz D, Förster DW, Wilting A. Phylogeography of red muntjacs reveals three distinct mitochondrial lineages. BMC Evol Biol. 2017;17:1–12. https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-017-0888-0.

    Article  Google Scholar 

  26. Pang JF, Kluetsch C, Zou XJ, Zhang AB, Luo LY, Angleby H, Ardalan A, Ekström C, Sköllermo A, Lundeberg J, Matsumura S. mtDNA data indicate a single origin for dogs south of Yangtze River, less than 16,300 years ago, from numerous wolves. Mol Biol Evol. 2009;26:2849–64. https://0-doi-org.brum.beds.ac.uk/10.1093/molbev/msp195.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Talukdar BK, Sharma A, Guleria H, Gupta M, Dudhwa’s rhinos—a plan for their growth & secured future. 2012. https://0-doi-org.brum.beds.ac.uk/10.13140/RG.2.2.10721.53601.

  28. Xu X, Janke A, Arnason U. The complete mitochondrial DNA sequence of the greater Indian rhinoceros, Rhinoceros unicornis, and the phylogenetic relationship carnivora, perissodactyla, and artiodactyla (+Cetacea). Mol Biol Evol. 1996;13:1167–73. https://0-doi-org.brum.beds.ac.uk/10.1093/oxfordjournals.molbev.a025681.

    Article  CAS  PubMed  Google Scholar 

  29. Steiner CC, Ryder OA. Molecular phylogeny and evolution of the Perissodactyla. Zool J Linn Soc. 2011;163:1289–303. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1096-3642.2011.00752.x.

    Article  Google Scholar 

  30. Vidya TNC. Evolutionary history and population genetic structure of Asian elephants in India. Indian J Hist Sci. 2016;51(2):391–405. https://0-doi-org.brum.beds.ac.uk/10.16943/ijhs/2016/v51i2.2/48453.

    Article  Google Scholar 

  31. Kumar A, Ghazi MGU, Hussain SA, Bhatt D, Gupta SK. Mitochondrial and nuclear DNA based genetic assessment indicated distinct variation and low genetic exchange among the three subspecies of swamp deer (Rucervus duvaucelii). Evol Biol. 2017;44(1):31–42. https://0-doi-org.brum.beds.ac.uk/10.1007/s11692-016-9387-2.

    Article  Google Scholar 

  32. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214. https://0-doi-org.brum.beds.ac.uk/10.1186/147121487214.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Ritchie AM, Lo N, Ho SYW. The impact of the tree prior on molecular dating of data sets containing a mixture of inter- and intraspecies sampling. Syst Biol. 2017;66:413–25. https://0-doi-org.brum.beds.ac.uk/10.1093/sysbio/syw095.

    Article  PubMed  Google Scholar 

  34. Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006;4: e88. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pbio.0040088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Yoder A, Yang Z. Estimation of primate speciation dates using local molecular clocks. Mol Biol Evol. 2000;17:1081–90. https://0-doi-org.brum.beds.ac.uk/10.1093/oxfordjournals.molbev.a026389.

    Article  CAS  PubMed  Google Scholar 

  36. Subramanian S, Denver DR, Millar CD, Heupink T, Aschrafi A, Emslie SD, Baroni C, Lambert DM. High mitogenomic evolutionary rates and time dependency. Trends Genet. 2009;2511:482–6. https://0-doi-org.brum.beds.ac.uk/10.1016/j.tig.2009.09.005.

    Article  CAS  Google Scholar 

  37. Knaus BJ, Cronn R, Liston A, Pilgrim K, Schwartz MK. Mitochondrial genome sequences illuminate maternal lineages of conservation concern in a rare carnivore. BMC Ecol. 2011;11(1):1–14. https://0-doi-org.brum.beds.ac.uk/10.1186/1472-6785-11-10.

    Article  CAS  Google Scholar 

  38. Ho SYW, Shapiro B. Skyline-plot methods for estimating demographic history from nucleotide sequences. Mol Ecol Resour. 2011;11(3):423–34. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1755-0998.2011.02988.x.

    Article  PubMed  Google Scholar 

  39. Moritz C. Defining ‘evolutionarily significant units’ for conservation. Trend Ecol Evol. 1994;9(10):373–5. https://0-doi-org.brum.beds.ac.uk/10.1016/0169-5347(94)90057-4.

    Article  CAS  Google Scholar 

  40. Crandall KA, Bininda-Emonds OR, Mace GM, Wayne RK. Considering evolutionary processes in conservation biology. Trend Ecol Evol. 2000;15(7):290–5. https://0-doi-org.brum.beds.ac.uk/10.1016/S0169-5347(00)01876-0.

    Article  CAS  Google Scholar 

  41. Kitchener AC, Breitenmoser-Würsten CH, Eizirik E, Gentry A, Werdelin L, Wilting A, Yamaguchi N, Abramov AV, Christiansen P, Driscoll C, Duckworth JW, Johnson W, Luo SJ, Meijaard E, O’Donoghue P, Sanderson J, Seymour K, Bruford M, Groves C, Hoffmann M, Nowell K, Timmons Z, Tobe S. A revised taxonomy of the Felidae. The final report of the Cat Classification Task Force of the IUCN/SSC Cat Specialist Group. Cat News Special Issue 2017;11:80.

  42. Karanth KK, Nichols JD, Karanth KU, Hines JE, Christensen NL Jr. The shrinking ark: patterns of large mammal extinctions in India. Proc R Soc B. 2010;277:1971–9. https://0-doi-org.brum.beds.ac.uk/10.1098/rspb.2010.0171.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Rookmaker K, Sharma A, Bose J, Thapa K, Dutta D, Jeffries B, Williams AC, Ghose D, Gupta M, Tornikoski S. The greater one-horned rhino: past, present and future. WWF, Gland, Switzerland, WWF. 2016.

  44. Talukdar BK, Bonal BS, Sharma A, Sarma K, Singh SP. Re-introduction of greater one-horned rhino in Manas National Park, Assam, India—under the Indian Rhino Vision 2020. Global re-introduction perspectives: 2013. further case studies from around the globe, 2020;164. https://www.researchgate.net/publication/346443488_Reintroduction_of_greater_one_horned_rhino_in_Manas_National_Park_Assam_India_-_under_Indian_Rhino_Vision_2020.

  45. Weeks AR, Dorian M, Thavornkanlapachai R, Taylor HR, White NE, Weiser EL, Heinze. Conserving and enhancing genetic diversity in translocation programs. In: Advances in reintroduction biology of Australian and New Zealand Fauna. Melbourne: CSIRO Publishing; 2017. p. 127–40.

    Google Scholar 

  46. Rookmaaker L. Former distribution of the Indian rhinoceros (Rhinoceros unicornis) in India and Pakistan. J Bombay Nat Hist Soc. 1983;80:55–563.

    Google Scholar 

  47. Talukdar BK, Emslie R, Bist SS, Choudhury A, Ellis S, Bonal BS, Malakar MC, Talukdar BN, Barua M. Rhinoceros unicornis: the IUCN red list of threatened species. 2008. http://www.iucnredlist.org/details/19496/0.

  48. Pant G, Maraseni T, Apan A, Allen BL. Trends and current state of research on greater one-horned rhinoceros (Rhinoceros unicornis): a systematic review of the literature over a period of 33 years (1985–2018). Sci Total Environ. 2020;710: 136349. https://0-doi-org.brum.beds.ac.uk/10.1016/j.scitotenv.2019.136349.

    Article  CAS  PubMed  Google Scholar 

  49. Hostetler JA, Onorato DP, Nichols JD, Johnson WE, Roelke ME, O’Brien SJ, Jansen D, Oli MK. Genetic introgression and the survival of Florida panther kittens. Biol Conserv. 2010;143:2789–96. https://0-doi-org.brum.beds.ac.uk/10.1016/j.biocon.2010.07.028.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Johnson WE, Onorato DP, Roelke ME, Land DE, Cunningham M, Belden C, McBride R, Jansen D, Lotz M, Shindle D, Howard J, Wildt DE, Penfold LM, Hostetler JA, Oli MK, O’Brien SJ. Genetic restoration of the Florida panther. Science. 2010;329:1641–5. https://0-doi-org.brum.beds.ac.uk/10.1126/science.1192891.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. O’Brien SJ, Johnson WE, Driscol CA, Dobrynin P, Marker L. Conservation genetics of the cheetah: lessons learned and new opportunities. J Hered. 2017. https://0-doi-org.brum.beds.ac.uk/10.1093/jhered/esx047.

    Article  PubMed  Google Scholar 

  52. Sharma SP, Ghazi MG, Katdare S, Dasgupta N, Mondol S, Gupta SK, Hussain SA. Microsatellite analysis reveals low genetic diversity in managed populations of the critically endangered gharial (Gavialis gangeticus) in India. Sci Rep. 2021;11:1–10. https://0-doi-org.brum.beds.ac.uk/10.1038/s41598-021-85201-w.

    Article  CAS  Google Scholar 

  53. Laurie A. Behavioural ecology of the greater one-horned rhinoceros (Rhinoceros unicornis). J Zool. 1982;196:307–41. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1469-7998.1982.tb03506.x.

    Article  Google Scholar 

  54. Johnsingh AJT, Manjrekar N. Mammals of South Asia, vol. II. India: Universities Press; 2016.

    Google Scholar 

  55. Ghosh T, Sharma A, Mondol S. Optimisation and application of a forensic microsatellite panel to combat greater-one horned rhinoceros (Rhinoceros unicornis) poaching in India. Forensic Sci Int Genet. 2021;52: 102472. https://0-doi-org.brum.beds.ac.uk/10.1016/j.fsigen.2021.102472.

    Article  CAS  PubMed  Google Scholar 

  56. Biswas S, Bhatt S, Paul S, Modi S, Ghosh T, Habib B, Nigam P, Talukdar G, Pandav B, Mondol S. A practive faeces collection protocol for multidisciplinary research in wildlife science. Curr Sci. 2019;116:11. https://0-doi-org.brum.beds.ac.uk/10.18520/cs/v116/i11/1878-1885.

    Article  CAS  Google Scholar 

  57. Hassanin A, Ropiquet A, Couloux A, Cruaud C. Evolution of the mitochondrial genome in mammals living at high altitude: new insights from a study of the tribe Caprini (Bovidae, Antilopinae). J Mol Evol. 2009;68:293–310. https://0-doi-org.brum.beds.ac.uk/10.1007/s00239-009-9208-7.

    Article  CAS  PubMed  Google Scholar 

  58. Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;13:1870–4. https://0-doi-org.brum.beds.ac.uk/10.1093/molbev/msw054.

    Article  CAS  Google Scholar 

  59. Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, Pütz J, Middendorf M, Stadler PF. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013;69:313–9. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ympev.2012.08.023.

    Article  PubMed  Google Scholar 

  60. Greiner S, Lehwark P, Bock R. Organellar genome DRAW (OGDRAW) version 1.3.1: expanded toolkit for the graphical visualization of organellar genomes. Nucleic Acids Res. 2019;47:W59–64. https://0-doi-org.brum.beds.ac.uk/10.1093/nar/gkz238.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Willerslev E, Gilbert MTP, Binladen J, Ho SYW, Campos PF, Ratan A, Tomsho LP, da Fonseca RR, Sher A, Kuznetsova TV, Nowak-Kemp M, Roth TL, Miller W, Schuster SC. Analysis of complete mitochondrial genomes from extinct and extant rhinoceroses reveals lack of phylogenetic resolution. BMC Evol Biol. 2009;9:95. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2148-9-95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Xu X, Arnason U. The complete mitochondrial DNA sequence of the white rhinoceros, Ceratotherium simum, and comparison with the mtDNA sequence of the Indian rhinoceros, Rhinoceros unicornis. Mol Phylogenet Evol. 1997;7:189–94. https://0-doi-org.brum.beds.ac.uk/10.1006/mpev.1996.0385.

    Article  CAS  PubMed  Google Scholar 

  63. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2. https://0-doi-org.brum.beds.ac.uk/10.1093/bioinformatics/btp187.

    Article  CAS  PubMed  Google Scholar 

  64. Bandelt H, Forster P, Rohl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48. https://0-doi-org.brum.beds.ac.uk/10.1093/oxfordjournals.molbev.a026036.

    Article  CAS  PubMed  Google Scholar 

  65. Leigh JW, Bryant D. PopART: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6:1110–6. https://0-doi-org.brum.beds.ac.uk/10.1111/2041-210X.12410.

    Article  Google Scholar 

  66. Corander J, Tang J. Bayesian analysis of population structure based on linked molecular information. Math Biosci. 2007;205:19–31. https://0-doi-org.brum.beds.ac.uk/10.1016/j.mbs.2006.09.015.

    Article  PubMed  Google Scholar 

  67. Excoffier L, Laval G, Schneider S. Arlequin ver. 3.0: an integrated software package for population genetics data analysis. Evol Bioinform. 2005;1:47–50.

    Article  CAS  Google Scholar 

  68. Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics. 2001;17:754–5. https://0-doi-org.brum.beds.ac.uk/10.1093/bioinformatics/17.8.754.

    Article  CAS  PubMed  Google Scholar 

  69. Harley EH, de Waal M, Murray S, O’Ryan C. Comparison of whole mitochondrial genome sequences of northern and southern white rhinoceroses (Ceratotherium simum): the conservation consequences of species definitions. Conserv Genet. 2016;17:1285–91. https://0-doi-org.brum.beds.ac.uk/10.1007/s10592-016-0861-2.

    Article  Google Scholar 

  70. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772–772. https://0-doi-org.brum.beds.ac.uk/10.1038/nmeth.2109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Eisenmann V. Origins, dispersals, and migrations of Equus (Mammalia, Perissodactyla). Cour Forschungsinst Senckenberg. 1992;153:161–70.

    Google Scholar 

  72. Orlando L, Ginolhac A, Zhang G, Froese D, Albrechtsen A, Stiller M, Schubert M, Cappellini E, Petersen B, Moltke I, Johnson PLF, Fumagalli M, Vilstrup JT, Raghavan M, Korneliussen T, Malaspinas A-S, Vogt J, Szklarczyk D, Kelstrup CD, Vinther J, Dolocan A, Stenderup J, Velazquez AMV, Cahill J, Rasmussen M, Wang X, Min GD, Zazula J, Seguin-Orlando A, Mortensen C, Magnussen K, Thompson JF, Weinstock J, Gregersen K, Røed KH, Eisenmann V, Rubin CJ, Miller DC, Antczak DF, Bertelsen MF, Brunak S, Al-Rasheid KAS, Ryder O, Andersson L, Mundy J, Krogh A, Gilbert MTP, Kjær K, Sicheritz-Ponten T, Jensen LJ, Olsen JV, Hofreiter M, Nielsen R, Shapiro B, Wang J, Willerslev E. Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse. Nature. 2013;499:74–8. https://0-doi-org.brum.beds.ac.uk/10.1038/nature12323.

    Article  CAS  PubMed  Google Scholar 

  73. Hooijer DA. Prehistoric and fossil rhinoceroses from the Malay Archipelago and India. Zool Meded. 1946;26:1–138.

    Google Scholar 

  74. Prothero DR, Schoch RM. The evolution of Perissodactyls (No. 15). Oxford University Press on Demand; 1989

  75. Benton MJ, Donoghue PCJ. Paleontological evidence to date the tree of life. Mol Biol Evol. 2007;24:26–53. https://0-doi-org.brum.beds.ac.uk/10.1093/molbev/msl150.

    Article  CAS  PubMed  Google Scholar 

  76. Rambaut A, Drummond AJ. Tracer. 2007. http://beast.bio.ed.ac.uk/tracer.

  77. Helfrich P, Rieb E, Abrami G, Lüking A, Mehler A. TREEANNOTATOR: versatile visual annotation of hierarchical text relations. In: Proceedings of the 11th edition of the language resources and evaluation conference, LREC 2018, Miyazaki, Japan.

  78. Rambaut A. FigTree v1. 3.1. 2009. http://tree.bio.ed.ac.uk/software/figtree/.

Download references

Acknowledgements

We thank the Ministry of Environment, Forest and Climate Change, Government of India for all support to implement this project. We thank Forest Departments of Assam, West Bengal and Uttar Pradesh for necessary permits and field support during sampling. We thank the WWF-India team for their field and logistic support. We appreciate the help from Mirza Ghazanfarullah Ghazi regarding marker aliquots during initial standardizations. Our thanks to Ankit Pacha, Bhim Singh, Surya P. Sharma and Dr. Kunal Arekar for their technical inputs in analysis. We acknowledge Dr. Jahnavi Joshi for her inputs in phylogenetic analysis. We thank the Director, Dean, Research Coordinator and Nodal Officer of the Wildlife Forensics and Conservation Genetics Cell for their support. We would also like to acknowledge three anonymous reviewers for their constructive comments.

Funding

This research was funded by MoEF&CC, Government of India and WWF India. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

Author information

Authors and Affiliations

  1. Animal Ecology and Conservation Biology Department, Wildlife Institute of India, Chandrabani, Dehradun, Uttarakhand, 248001, India

    Tista Ghosh, Shrewshree Kumar, Kirtika Sharma & Samrat Mondol

  2. World Wide Fund for Nature-India, 172B Lodhi Estate, New Delhi, 110003, India

    Parikshit Kakati & Amit Sharma

Authors
  1. Tista Ghosh
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Shrewshree Kumar
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Kirtika Sharma
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Parikshit Kakati
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Amit Sharma
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Samrat Mondol
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

TG and SM designed the study. AS and SM raised necessary funds and acquisition of permissions to conduct this study. Sample collection was done by TG, SK and PK. TG, SK and KS were involved in data curation and generation, whereas analysis and data visualisation were done by TG and SK. TG and SM wrote the manuscript with inputs from SK, AS, PK and KS. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Tista Ghosh or Samrat Mondol.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Details of the primers used in sequencing Indian rhino whole mitogenome and phylogeography data. Table S2. Sample (both tissue and dung) of Indian rhinos used in this study. A total of 111 individual rhino samples (72 tissue and 39 dung samples, respectively) were used in this study. Out of the 72 tissue samples, 16 samples representing different parks were used to generate the whole mitogenome data. Table S3. Mitogenome organization in Rhinoceros unicornis. Codons respective to each tRNA are mentioned in parenthesis. Table S4. Comparative analysis results of genetic diversity indices among five extant rhino species. Table S5. Details of the variable sites based on concatenated sequence of 2531bp of Indian one horned rhino mtDNA haplotypes.

Additional file 2: Figure S1.

Whole mitogenome organisation and annotation of Rhinoceros unicornis. Figure S2. Gene-wise comparative analyses of polymorphism indices (S, Hd and π) for all extant rhino species. Figure S3. Estimation of mitogenome mutation rate for Indian rhino using Caballine (Zebra, Donkey, Asiatic Wild Ass and Horse) as outgroup. Internal node calibration points are presented in italics (bold). Estimated node ages and branch mutation rates with available reference in literature is marked with * [1, 8, 10, 29].

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ghosh, T., Kumar, S., Sharma, K. et al. Consideration of genetic variation and evolutionary history in future conservation of Indian one-horned rhinoceros (Rhinoceros unicornis). BMC Ecol Evo 22, 92 (2022). https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-022-02045-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-022-02045-2

Keywords

BMC Ecology and Evolution

ISSN: 2730-7182

Contact us