Skip to main content
Download PDF

Geography vs. past climate: the drivers of population genetic structure of the Himalayan langur

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

A Correction to this article was published on 23 August 2022

This article has been updated

Abstract

Background

Contemporary species distribution, genetic diversity and evolutionary history in many taxa are shaped by both historical and current climate as well as topography. The Himalayas show a huge variation in topography and climatic conditions across its entire range, and have experienced major climatic fluctuations in the past. However, very little is known regarding how this heterogenous landscape has moulded the distribution of Himalayan fauna. A recent study examined the effect of these historical events on the genetic diversity of the Himalayan langurs in Nepal Himalaya. However, this study did not include the samples from the Indian Himalayan region (IHR). Therefore, here we revisit the questions addressed in the previous study with a near complete sampling from the IHR, along with the samples from the Nepal Himalaya. We used the mitochondrial Cytochrome-b (Cyt-b, 746 bp) region combined with multiple phylogeographic analyses and palaeodistribution modelling.

Results

Our dataset contained 144 sequences from the IHR as well as the Nepal Himalaya. Phylogenetic analysis showed a low divergent western clade nested within high divergent group of eastern lineages and in the network analysis we identified 22 haplotypes over the entire distribution range of the Himalayan langurs. Samples from the Nepal Himalaya showed geographically structured haplotypes corresponding to different river barriers, whereas samples from IHR showed star-like topology with no structure. Our statistical phylogeography analysis using diyABC supported the model of east to west colonisation of these langurs with founder event during colonisation. Analysis of demographic history showed that the effective population size of the Himalayan langurs decreased at the onset of last glacial maximum (LGM) and started increasing post LGM. The palaeodistribution modelling showed that the extent of suitable habitat shifted from low elevation central Nepal, and adjoining parts of north India, during LGM to the western Himalaya at present.

Conclusion

The current genetic diversity and distribution of Himalayan langurs in the Nepal Himalaya has been shaped by river barriers, whereas the rivers in the IHR had relatively less time to act as a strong genetic barrier after the recent colonisation event. Further, the post LGM expansion could have had confounding effect on Himalayan langur population structure in both Nepal Himalaya and IHR.

Peer Review reports

Background

Accumulation of genetic variation across different populations of a species could be attributed to changes in geography and climate via a combined effect of genetic drift, gene flow and selection [1,2,3]. Geographical barriers such as mountain ranges and rivers; environmental barriers; and certain behavioural aspects of the taxa such as philopatry, dispersal ability, etc. can prevent gene flow, and over time create population genetic structures [4, 5].

Rivers are considered as long-term barriers to gene flow and therefore plays an important role in speciation and divergence [6,7,8]. Wallace [9] noted that many south American primate species are separated by large rivers in the Amazon basin, and that different species of primates are restricted to opposite banks of the river and never cross it. The ability of a river to act as a barrier depends on various physical attributes of the river such as the size, the amount of discharge and speed of the water flow; further the dispersal ability, ecology and body size of the organism also plays an important role [10, 11]. The effects of rivers as barriers have been studied in many taxa of plants, amphibians, reptiles, mammals and birds [12,13,14,15,16,17,18,19,20,21]. In primates, these studies of riverine barrier effect have mostly been confined to Amazonian and African taxa [4, 10, 11, 22,23,24], with only a handful of studies in Asian primates [25,26,27,28,29]. Apart from the geographical barriers, Climatic fluctuations are also known to have an impact on the population genetic structure and demography [30,31,32,33]. Quaternary glaciation periods, such as the last glacial maximum (LGM), substantially influenced how genetic variation is spatially distributed in many species across the globe [31, 34]. Repeated periods of glaciation, especially in higher latitudes, resulted in range expansion—in the case of species adapted to cold climate; and range contraction—in the species adapted to warm climate [31, 34,35,36,37,38,39,40,41,42,43,44,45]. In the latter case, species may undergo range expansion and exponential population growth after the end of the glacial period [46, 47].

The Himalayan range is a globally recognised biodiversity hotspot with high levels of endemism and a complex geography, topography and climate [48]. It extends 2500 kms in length from west to east and its width varies between 350 km in the west and 150 km in the east [49]. Given these, there are very few studies testing the effect of such a complex topography and climatic fluctuations on the distribution and demography of species found in the Himalayas [28, 29, 50, 51]. Himalayan langur, (Primate: Colobinae) is a widely distributed alloprimate found in most parts of the Himalayas. It is distributed in Nepal and parts of India, Pakistan and Bhutan [52], the Black mountain range and the Sunkosh river in Bhutan forms the easternmost distribution limit for this species [53]. Until recently, there was much ambiguity in the taxonomy of these langurs owing to presence of multiple classification schemes, however IUCN recognised three species—Semnopithecus ajax, Semnopithecus schistaceus and Semnopithecus hector, as per the classification scheme proposed by Groves [54]. Nevertheless, there was no consensus within the scientific community on the use of any one classification scheme and the names proposed by all the schemes were considered to be valid and were used by different studies [28, 55,56,57,58,59,60]. This use of different taxonomic names for the same species can create serious issues while interpreting results from these studies, therefore, it was important to resolve the taxonomy of this group. In this context, a recent systematic study resolved the taxonomy of these langurs using an integrative approach based on three lines of evidence from molecular, morphological and ecological data. The authors have now subsumed the three recognised species into a single widespread species Semnopithecus schistaceus Hodgson, 1840 [52].

Another study investigated the role of Himalayan rivers as well as past climate in shaping the distribution of Himalayan langur in Nepal Himalaya [28]. They found that Himalayan langur populations (they used the name S. entellus in their paper) in the Nepal Himalayas were isolated into six major clades by the snow-fed Himalayan rivers suggesting the role of these rivers as a barrier to gene flow. Further, their demographic analysis showed that the Himalayan langur populations in Nepal started declining with the onset of the last glacial maximum (LGM) i.e., ~ 22,000 years before present (YBP). For their study, they used two mitochondrial genes—Cytochrome b (Cyt b) and hyper variable region (HVR) I. However, this study did not include any samples from the Indian Himalayan region (IHR) which constitutes a major part of the Himalayan langur distribution zone.

In this study, we revisit the questions addressed by Khanal et al. [28] with near complete sampling from India. Along with the samples from the Nepal Himalaya [28], we have expanded the sampling into the Indian Himalayan region and included three river valleys from this region—Sutlej river valley, Bhagirathi river valley, and Mahakali/Sharda (hereafter Mahakali) river valley. Given this extensive sampling, covering the entire distribution range of the Himalayan langur in this study, we first investigate the role of major Himalayan rivers in shaping the population structure of Himalayan langurs. For this, we have included all the river valleys studied by Khanal et al. [28] (except the river Trishuli since it did not act as barrier for gene flow as per [28]) along with three river valleys from the IHR west of Nepal. Additionally, given the wider sampling of Himalayan langurs in this study we revisit the role of past climate, specifically the last glacial maximum (LGM), in shaping the contemporary distribution of the Himalayan langur. Finally, we investigate the purported westward expansion of the Himalayan langur from Nepal and Sikkim. In the molecular phylogenetic tree of Arekar et al. [52], a low divergent western clade was nested within high divergent eastern lineages suggesting a potentially east to west dispersal. Therefore, we considered the western clade as one population namely the western metapopulation (WMP) and the high divergent eastern lineages as another population i.e., the eastern metapopulation (EMP). To test the east to west dispersal hypothesis, we used a model-based hypothesis testing in a Bayesian framework using Approximate Bayesian computation (ABC) analysis implemented in DIY ABC v2.1.0 [61].

Results

Phylogenetic analysis

The final dataset contained 746 bp of 145 Himalayan langur sequences, out of these 51 were sequenced in this study; 25 sequences were downloaded from Arekar et al. [52]; 67 sequences from Khanal et al. [28]; 1 sequence from Karanth et al. [62] and one S. entellus sequence was used as an outgroup (Additional file 4: Table S1) [28, 52, 62]. HKY + G substitution model was selected for this analysis. In both the Bayesian (Fig. 1) and ML tree (Additional file 1: Fig. S1), we found a well-supported clade—the western clade, consisting of haplotypes from the western Himalayas (west of the river KaliGandaki), which was nested within samples from eastern Himalayas (east of the river KaliGandaki) i.e., the eastern lineages. Both the trees showed similar topology wherein all the major clades were retrieved. The only exception was in the Bayesian tree (Fig. 1), the MH271128 and MH271129 sequences were placed within the clade containing samples from the Indian Himalayan region (IHR), whereas in the ML tree (Additional file 1), these two lineages were sister to the IHR clade. The position of the MH271121 clade was also different in the Bayesian and ML trees. Our divergence dating analysis showed that the western clade (i.e. WMP) diverged from the eastern lineages (i.e. EMP) at around 0.52 mya (CI 0.28–0.78). Further, the population in IHR diverged from the rest of the Himalayan langurs at 0.29 mya (CI 0.16–0.44) (Additional file 2).

Fig. 1
figure 1

Bayesian phylogeny of Himalayan langur for mitochondrial cytochrome b (Cyt-b) gene. Numbers at the node indicate Bayesian posterior probability (BPP) values. Node support values > 0.75 are shown. Node support values are shown only for major clades. The colours for the tip taxa correspond to the colours of the 10 populations in the haplotype network in Fig. 2

Full size image

Phylogeographic analyses

For these analyses, we removed the outgroup sequence from our dataset and the sequence length was shortened to 645 bp due to removal of sites with missing data. So, the final dataset for these analyses contained 144 sequences.

Network analysis, genetic diversity, and population genetic structure

The network analysis of 144 sequences yielded 22 haplotypes across the distribution zone of the Himalayan langurs (Fig. 2). Two clusters can be observed in the network, one consisting of the WMP (H1–H13) i.e., all the sequences from the western Himalayan region (west of the river KaliGandaki), and the other consisting of the EMP (H14–H22) i.e., all the sequences from the eastern Himalayan region (east of the river KaliGandaki). These two clusters are spatially delineated by the river Kali Gandaki. Within the western population, one high frequency haplotype (H7, n = 29) was observed which consisted majority of the individuals from the western region. The haplotypes H10 and H11 were separated from the rest of the samples of the IHR by river Mahakali; (these two haplotypes correspond to the WB clade as labelled in Khanal et al. [28]). The haplotypes H12 and H13 were also seen forming a separate cluster separated from the rest by river Karnali, these two haplotypes correspond to the clade WA in Khanal et al. Other than these, the samples from the IHR showed very little structure across the river valleys. On the contrary, the eastern population showed highly structured haplotypes separated by the river valleys in the Nepal Himalayas as shown in Fig. 3. The river Marsyagandi separated the haplotype H18 (which corresponds to clade CC in Khanal et al. [28]) from the haplotypes H19 and H20 (these haplotypes correspond to clade CB in Khanal et al. [28]). And the river Budhi Gandaki separated the haplotypes H19 and H20 from the haplotypes H21 and H22 (these haplotypes correspond to clade CA in Khanal et al. [28]). Further, the rivers Arun and Tamor separated the haplotypes H16 and H17, respectively. Here the haplotype H16 corresponds to the clade EA in Khanal et al. [28]. In this study, the haplotype H14 is placed in the population Marsyagandi_BudhiGandaki (My–Bg) based on the geographical proximity however, in the haplotype network it can be seen as an outlier; and it corresponds to the haplotype H19 in Khanal et al. [28]. This pattern seen in the haplotype network analysis is similar to the results from our phylogenetic analysis (Fig. 1).

Fig. 2
figure 2

Haplotype network for the Himalayan langurs using Cyt-b sequence. The coloured circles represent extant haplotypes while the black circles are inferred intermediate haplotypes not sampled in this study. The sizes of the haplotypes are proportional to their frequency i.e., the number of individuals that constitute those haplotypes. The bars on the link between the circles stand for the number of substitutions between those haplotypes. The dotted red and green lines indicate the two larger groups—WMP and EMP, respectively

Full size image
Fig. 3
figure 3

Diagram showing the 10 populations of Himalayan langurs separated by the nine river valleys; 1—River Sutlej, 2—River Bhagirathi, 3—River Mahakali, 4—River Karnali, 5—River KaliGandaki, 6—River Marsyagandi, 7—River BudhiGandaki, 8—River Arun, 9—River Tamor. The colours of these populations correspond with the haplotype colours in Fig. 2. The yellow borderline demarcates the country of Nepal, its extent is also shown in the map using a square bracket

Full size image

The nucleotide diversity (π) of the My–Bg population was highest (0.015), in fact it was equivalent to all the populations combined together; the haplotype diversity (H) was also high (0.8) (Table 1). The populations to the west and to the east of the My–Bg has low π and H values. The genetic diversity data for three populations; KaliGandaki_Marsyagandi (Kg–My), Arun_Tamor (A–T) and EastTamor (E-T) was not obtained due to the absence of polymorphic sites in the nucleotide sequences.

Table 1 DNA polymorphism and genetic diversity of different populations of the Himalayan langurs
Full size table

The pairwise FST values between populations across WMP and EMP were higher than between populations within WMP or EMP (Table 2). Further, the FST values were low among populations that were geographically closer whereas they were higher in spatially distant populations.

Table 2 Comparisons of pairs of population samples—population pairwise Fst (below diagonal), average number of pairwise difference between populations (above diagonal) and average number of pairwise differences within populations (diagonal elements in bold) among the different populations of Himalayan langurs calculated by distance method
Full size table

Demographic history

Our mismatch distribution analyses for populations in the WMP and the EMP, each produced different patterns (Fig. 4). The EMP yielded a multimodal distribution of the observed data (Fig. 4I). Within the EMP, the population My–Bg showed a multimodal distribution (Fig. 4F); all the other population either showed a unimodal distribution, or a bimodal distribution closer to y-axis.

Fig. 4
figure 4

Results of mismatch distribution carried out on Cyt-b data obtained from the Himalayan langurs. The graphs show pairwise differences between sequences (X-axes) plotted against the frequency of those differences to generate the distributions. The blue dotted line indicates expected distribution under a population growth-decline model. WMP western metapopulation, EMP eastern metapopulation. The panels A–I represent individual mismatch graphs for different populations, the population names are mentioned on the top right corner of each graph

Full size image

The test of neutrality to check for excess of rare mutations as evidence for population expansion was not significant, except for the population BudhiGandaki_Arun (Bg–A) which showed a significantly negative Tajima’s D value (− 1.844, p < 0.05) (Table 1). The Bayesian skyline plot analysis for the overall dataset showed that the effective population size (Ne) of the Himalayan langurs decreased during the LGM, and it started increasing post LGM ultimately reaching its current population size (Fig. 5). The separate Bayesian skyline plots for WMP and EMP were not informative.

Fig. 5
figure 5

Bayesian skyline plot reconstructed using Cyt-b gene for the Himalayan langur (S. schistaceus). This plot shows the changes in effective population size (Ne) through time. X-axis is time in years before present and Y-axis is the estimated effective population size. The solid line indicates the median Ne; the grey shaded area shows the 95% highest posterior density (HPD) intervals; vertical red dotted lines denotes the range of LGM from 26.5 to 19 ka as per [63]

Full size image

Statistical phylogeography

For the ABC analysis we tested three scenarios (Fig. 6), Scenario 1 hypothesised that WMP originates from EMP with no change in the effective population size; Scenario 2 hypothesised a founder event in WMP; and Scenario 3 assumed bottleneck in the ancestral population of WMP and EMP. Further details are explained in “Methods” section. Here, the second scenario was selected as the scenario of choice with highest posterior probabilities in all the summary statistics used. In this scenario, we hypothesised a founder event where few individuals from the EMP colonised the western Himalayas and the population size eventually increased to the current state (Fig. 6, scenario 2). Model checking too supported scenario 2, since among the three scenarios the simulated dataset under scenario 2 were closest to the observed values.

Fig. 6
figure 6

Alternative scenarios of demographic history of the Himalayan langur (S. schistaceus) tested using ABC analysis in DIYABC. Detailed explanation is mentioned in “Demographic history” in “Methods” section. Time has been measured in generations before present. The table alongside this figure contains details of the model specified, prior distributions for the demographic parameters and the mutation model parameters. WMP western metapopulation, EMP eastern metapopulation

Full size image

Niche modelling using past climate layers

The AUC values for the training and test data for the Himalayan langur dataset were 0.9772 and 0.9702, respectively, indicating that the potential distribution of the Himalayan langurs fits well with our data. Precipitation of the driest quarter (Bio 17) had the highest contribution to the model (45.1%) followed by annual mean temperature (Bio1; 26.2%) and precipitation seasonality (Bio15; 8.4%). The response curves reveal that Bio17 value of above 200, Bio 1 value between 80 and 120 and Bio 15 value in the range of 115–121 seem to be an ideal habitat for Himalayan langurs (Additional file 3). According to the palaeodistribution model the distribution of Himalayan langurs in the LGM was more towards south especially in the lowland Terai region of central Nepal and adjoining parts of northern India; however, the probability of distribution was moderate. Further, the probability of distribution of these langurs in the western Himalaya was higher for current time than during LGM (Fig. 7).

Fig. 7
figure 7

Ecological niche modelling projections of the Himalayan langur (Semnopithecus schistaceus). Top panel shows current distribution and the bottom panel shows potential distribution during LGM

Full size image

Discussion

After decades of taxonomic confusion, the Himalayan langurs were recently assigned to a single species Semnopethicus schistaceus [52]. These langurs have a wide distribution in the Himalayas and therefore are a perfect model system to study the effect of landscape heterogeneity on genetic diversity and population genetic structure. The Himalayan langurs are threatened by habitat destruction—a result of human mediated climate change and human population expansion [64]. A recent study, using species distribution modelling, predicted range contraction of ~ 64% by 2050 under the GCMs RCP 4.5 and 8.5 [65]. Here, we investigated the drivers of population demography and genetic structure of the Himalayan langurs. Results from this study can inform the future management and conservation actions for this species.

Genetic diversity and population structure

The study found that all the Himalayan langur samples, when pooled together, have a high haplotype diversity (Hd = 0.915) when compared to another temperate colobine Rhinopithecus brelichi (Hd = 0.457), and a low nucleotide diversity (π = 0.015) which is comparable to R. brelichi (π = 0.014) [66]. This combination of high haplotype diversity and low nucleotide diversity could be the result of population expansion after a recent bottleneck event or population expansion after a recent colonisation event [30]. To test for the riverine barrier hypothesis, we divided the Himalayan langur individuals into different populations demarcated by river valleys, here we found that the population My–Bg showed a high haplotype diversity (Hd = 0.8) and its nucleotide diversity was highest among all other populations (π = 0.015). This suggests that the My–Bg population which is demarcated by the river Marsyagandi and the river Budhi Gandaki in central Nepal, could either be a long-term stable population or it could also suggest that this population contains distinct haplotypes, probably originating from different sources. Rest of the populations showed a low π and Hd values, indicative of population expansion either after a recent bottleneck event or a recent colonisation event. Results from our pairwise FST analysis indicates that the WMP and EMP have been separated from each other longer than populations within WMP and EMP.

Do himalayan rivers act as barriers?

Our phylogenetic tree with wider sampling is very similar to the tree in Arekar et al. [52] wherein a western clade with low divergence (corresponding to WMP) was nested within a larger phylogeny consisting of highly divergent sequences from eastern region (corresponding to EMP). The network analysis (Fig. 2), generated two major clusters separated by the river Kali Gandaki—the western cluster and the eastern cluster which corresponds to WMP and EMP, respectively (as seen in Fig. 1). The western cluster showed a star-like phylogeny which suggests shallow genetic structure and recent demographic expansion [67]. Whereas, the eastern cluster showed higher divergence among different haplotypes. The rivers in the IHR did not cause any structuring among different populations, in that haplotypes are shared across rivers. However, the rivers in the Nepal Himalaya i.e., Mahakali, Karnali, Kali Gandaki, Marsyagandi, Budhi Gandaki, Arun and Tamor, shaped the high genetic structure among different populations. Two populations within the WMP falls in the Nepal Himalaya i.e. Mahakali_Karnali (M–Kr) and Karnali–KaliGandaki (Kr–Kg), these two populations shows structured haplotypes across the river Karnali. Given that WMP are younger than the populations in the east i.e., EMP, it is plausible that the rivers in the IHR did not have enough time to act as barriers. Rivers have been shown to act as barriers to gene flow across different taxa [22, 31, 68], however the ability of river to act as barrier depends on various factors such as its size, speed of the flow of water, body size of the organism, and the dispersal ability of the organism [22]. The Kali Gandaki river is characterised by its strong water currents and it forms the deepest gorge in the world, therefore not surprising it forms a barrier against geneflow between the two populations—WMP and EMP, of Himalayan langurs.

Effects of past glaciation events on demographic history

The effects of past climatic events, such as Pleistocene glaciation, on the distribution and demography of a species includes range contraction, range fragmentation, local extinction and resultant bottleneck [69, 70]. After the glaciation ends, species may undergo range expansion and exponential population growth [46, 47]. On the contrary, the habitats of taxa adapted to cold and dry conditions expands/persists during glaciation events [35,36,37, 44, 45].

In our analysis, the Tajima’s D, Fu’s F, and Ramos-Onsins and Rozas’s R2 statistic were not informative since the values for these indices were not significant in any of the population, except for the population Bg–A. However, it should be noted that these statistics by itself may not be very reliable for investigating the demographic history of a population and may depend on additional factors. So, for reliable predictions of demographic history, a combination of different indicators should be used [71]. In our mismatch distribution analysis, the EMP shows multimodal graph (Fig. 4I) which could either suggest that it is a long-term stable population or this population consists of haplotypes originating from more than one source. It should be noted that the mismatch distribution pattern is significantly affected by population structure [72, 73] and since haplotypes within EMP are highly structured (Fig. 2), we think that the mismatch analysis result here indicates the later. WMP also shows a multimodal graph but most likely it might be a result of population structure within WMP (Fig. 4H). We also performed mismatch analysis for individual populations (Fig. 4A–F); the three populations in the WMP which are distributed in the IHR showed a unimodal distribution which indicates population expansion (Fig. 4A–C). The populations M–Kr and Kr–Kg show a bimodal pattern which generally is an indicator of population under bottleneck, the data here is insufficient to support this result. And for the individual populations within the EMP, the populations My–Bg and Bg–A (Fig. 4F, G) showed a multimodal graph indicating that either these are long-term stable populations or the haplotypes they contain have originated from multiple source populations. Other populations within EMP did not show signatures of long-term stable population.

The Bayesian skyline plot indicated that the Himalayan langur population size was constant for a long period of time and started decreasing about ~ 25,000 years before present (YBP) i.e., after the onset of the last glacial maximum (LGM). After the LGM ended, the effective population size started increasing and reached the present-day estimate (Fig. 5). Our ABC analysis also supported the demographic expansion scenario in the WMP after it diverged from the EMP i.e., Scenario 2 (Fig. 6), which was supported by higher posterior probability. This suggests a founder effect during colonisation where a few individuals dispersed from Nepal Himalaya into IHR and later expanded to the current population size. However, our divergence dating analysis suggest that this demographic expansion event occurred before LGM.

Palaeodistribution modelling results suggested that during LGM, the Himalayan langur distribution was spread to lower elevation central Nepal and adjoining parts of India, although the probability of distribution was less. Precipitation of the driest quarter, annual mean temperature, and precipitation seasonality were the main contributing factors in defining the suitable habitat for the Himalayan langurs. Given the cold and dry conditions during LGM at high altitudes in the Himalayas and in contrast, warmer conditions with high precipitation [74, 75] at the lower elevations (a combined effect of south west monsoon and mid-latitude westerlies), suggests that the low elevation central Nepal and adjoining parts of India could have potentially acted as refugia for these langurs. After LGM ended, the Himalayan langur distribution moved northwards; it also facilitated the movement of these langurs from central Nepal into western Himalaya where a high probability for the distribution of these langurs can be seen in the current SDM (Fig. 7). However, in our SDM for the LGM layers, we can see that there were small pockets of high probability of distribution in the western region. The asynchronous glacial advances during the LGM in the Himalayas [75, 76], could explain these pockets of high probability of distribution in the western region, however, the phylogenetic and phylogeographic analysis of the Himalayan langurs do not show signatures of refugia in the western region.

Conclusion

We used multiple complementary methods to assess the genetic diversity, population structure and demographic history of the Himalayan langurs. In the network analysis, two distinct population clusters within the Himalayan langur were retrieved—WMP and EMP corresponding to the western clade and eastern lineages, respectively as described by Arekar et al. [52]. Interestingly, Himalayan river valleys in the IHR does not appear to affect the population subdivisions and distribution of genetic variation in the WMP, except for the populations M–Kr and Kr–Kg which showed structure across the river Karnali. The rivers in the Nepal Himalaya do act as a barrier to geneflow in the EMP. Our phylogenetic and statistical phylogeography analysis suggest a recent east to west dispersal of the Himalayan langurs. Given, the WMP is younger than EMP, it is likely that not enough time has elapsed for rivers in IHR to shape their genetic structure. Furthermore, post LGM population expansion of WMP might have also confounded their population structure. In this study, we have used only one mitochondrial Cyt-b gene; there might be a possibility that this gene did not capture enough variation, especially within the western population given the recent colonisation. Therefore, using markers with high substitution rates such as microsatellites and SNPs, could help us better understand the role of these river valleys in population subdivision within the Himalayan langurs.

Methods

Sample collection

We collected 176 fecal samples from 46 different locations in the Himalayas covering the Union Territory of Jammu and Kashmir (J&K), and the states of Himachal Pradesh (HP), Uttarakhand and Sikkim (Table S1 in [52]). Multiple samples were collected from each location. We successfully amplified a 775 bp (amplicon length) of mitochondrial Cytochrome b (Cyt-b) gene using the primer pairs Cytb_278F and Cytb_1052R. For details on DNA extraction, PCR amplification and sequencing see Arekar et al. [52].

Phylogenetic analysis and divergence dating

The sequence files obtained were viewed and edited manually in ChromasLite v2.01 (Technelysium Pty. Ltd.). Sequences generated in this study were combined with those generated by Khanal et al. [28] and aligned using MUSCLE algorithm [77] incorporated in MEGA v7 [78]. We used jModelTest 2.1.3 [79] to pick the best model of sequence evolution. Phylogenetic reconstruction was performed using Maximum Likelihood (ML) and Bayesian methods. ML analysis was performed in RAxML7.4.2 incorporated in raxmlGUI v1.3 [80]. 1000 replicates were performed to assess support for different nodes. We used MrBayes 3.2.2 [81] to perform the Bayesian analysis. Two parallel runs, with four chains each, were run for 5 million generations with sampling frequency set to every 1000 generations. Convergence between the two runs was determined based on standard deviation of split frequencies. The program Tracer v1.6 [82] was used to determine stationarity, an effective sample size (ESS) value of > 200 for each parameter was used as a cut-off for run length. The first 25% of trees were discarded as burn-in.

For divergence time estimation we used BEAST v2.6 [83]. The sequences used for this analysis and their accession numbers are shown in Additional file 4: Table S2. The input file was prepared using BEAUti v2.6 [83]. The data was partitioned into three codon-based partitions—codon position 1 with site model K80 + G, codon position 2 with site model HKY + I and codon position 3 with site model TN93 + I. The best substitution scheme and the model of sequence evolution was selected using PartitionFinder 1.1.0 [84]. Clock model was set to uncorrelated relaxed clock lognormal. For setting the clock rate parameter, we used the rate of substitution as estimated for Human mtDNA protein coding genes [85]. The first and second codon position rate was set to 8.8 × 10−9 substitutions per nucleotide per year and for the third codon position, rate was set to 1.9 × 10−8 substitutions per nucleotide per year. The analysis was run for 100 million generations with sampling frequency of 10,000.

Phylogeographic analyses

Network analysis, genetic diversity and population genetic structure

To explore the role of rivers in shaping the population genetic structure of Himalayan langurs sampling locations were divided into 10 different populations, corresponding to 10 regions in the Himalayas, demarcated by different river valleys—WestSutlej (W-S): sample locations west of the river Sutlej; Sutlej_Bhagirathi (S–B): sample locations between the rivers Sutlej and Bhagirathi; Bhagirathi_Mahakali (B–M): sample locations between the rivers Bhagirathi and Mahakali; Mahakali_Karnali (M–Kr): sample locations between the rivers Mahakali and Karnali; Karnali_KaliGandaki (Kr–Kg): sample locations between the rivers Karnali and Kali Gandaki; KaliGandaki_Marsyagandi (Kg–My): sample locations between the rivers Kali Gandaki and Marsyagandi; Marsyagandi_BudhiGandaki (My–Bg): sample locations between the rivers Marsyagandi and Budhi Gandaki; BudhiGandaki_Arun (Bg–A): sample locations between the rivers Budhi Gandaki and Arun; Arun_Tamor (A–T): sample locations between the rivers Arun and Tamor; EastTamor (E-T): sample locations east of the river Tamor (Fig. 3).

Genetic diversity indices including number of polymorphic sites (s), haplotype number (H), haplotype diversity (Hd) and nucleotide diversity (π) were calculated using DnaSP v6 [86] for each of these populations. Further, we used the median-joining (MJ) network [87] incorporated in PopART http://popart.otago.ac.nz/ [88] to build a haplotype network which graphically represents the relationship of each sample from different geographical locations. Population pairwise FST was calculated using Arlequin v3.5.2.2 and the statistical significance was tested by performing 10,000 permutations.

Demographic history

To infer the demographic history, we used multiple, complimentary methods. We calculated three summary statistics: Tajima’s D, Fu’s F, and Ramos-Onsins and Rozas’s R2 statistic for all the populations to understand the evolutionary history under different demographic scenarios using DnaSp v6 [86]. We also performed mismatch distribution analysis using the population growth-decline model to estimate the trends of population growth using DnaSp v6 [86]. We used the Bayesian Skyline Plot (BSP) analysis to estimate population size changes through time using BEAST v2.6 [83]. The best substitution model was selected using the BIC criterion in Modeltest v2.1.3 [79]. The clock model was set to strict clock and the mutation rate was set to 0.0178 mutations per site per million years [70]. The tree prior was set to BirthDeath Skyline Contemporary; and the Origin_BDSKY_Contemp.t prior was set as lognormal with mean of 2.7 and SD of 0.5. This prior sets the date of the root node, in this case the divergence between S. entellus and S. schistaceus. The analysis was run for 750 million generations with sampling every 1000 generations. We used Tracer v1.6 to check for stationarity by ensuring that the effective sample size (ESS) for each parameter was > 200. The BDSKY plot was visualised using a R script [89].

Statistical phylogeography

To understand the westward expansion of these langurs, we carried out a model based hypothesis testing in a Bayesian framework using Approximate Bayesian computation (ABC) analysis implemented in DIY ABC v2.1.0 [61]. ABC approximates posterior probabilities which it then uses to rank the different scenarios being considered. It first creates a prior distribution of parameter values by simulating large number of datasets under each scenario and then it uses a logistic regression method to estimate the posterior probability by picking scenarios, which are from among the simulated datasets, that are closest to the observed data [90]. For this analysis we tested three scenarios (Fig. 6), Scenario 1 hypothesised that WMP originates from EMP with no change in the effective population size; In Scenario 2 we hypothesise a founder event where few individuals from EMP colonised the western region, and the population eventually increased to the current size. and Scenario 3 assumed bottleneck in the ancestral population of WMP and EMP with increase in effective population size immediately after the divergence of the two populations. In the first two scenarios, EMP is considered as ancestral population because the phylogeny (Fig. 1) shows that WMP is nested within the EMP. The prior settings for the demographic model and the mutation model are shown in the table alongside Fig. 6. We estimated four one sample and four two sampled summary statistics. One million datasets were simulated for each scenario.

Niche modelling using past climate layers

For our second hypothesis, we implemented the Ecological niche modelling (ENM) approach. Here we used 217 occurrence records of the Himalayan langurs. Out of these, 104 records were from the field surveys conducted for this study, 58 occurrence records were obtained from previous studies [28, 91, 92] and 55 occurrence records were downloaded from GBIF (Global Biodiversity Information Facility) database (www.gbif.org). The occurrence records will be made available upon request. We used the MaxEnt algorithm [93] and 19 bioclimatic variables (www.worldclim.org) for the current (~ 1960–2000) and last glacial maximum (LGM) (~ 22,000 years before present, YBP) layers. Spatial resolution of LGM was resampled to 30 arcsec to match the current layers. These bioclimatic variables were clipped to the region from 68 °E to 97.4 °E and from 6.7 °N to 37 °N using ArcGIS 10.2.1. These clipped layers were then exported to ASCII format using QGIS 2.18.12. The 19 bioclimatic layers were tested for multicollinearity by calculating Pearson’s correlation coefficient (r). The layers with r ≤ |0.8| were selected for further analysis. Performance of Maxent depends on the choice of features and regularisation multiplier (RM) [93]. We tested 48 models, for the Himalayan langur dataset by employing different combinations of features and RM values (Additional file 3: Table S3) in MaxEnt v3.4.1 [93]. The model (a combination of features and RM) with the highest AUC value was selected as the best model.

ENM analysis was performed in Maxent v3.4.1 with the following modifications; Random test percentage was set to 30%, maximum number of background points was set to 10,000 and the replicates were set to 10 with replicated run type changed to Subsample. 5000 iterations were performed with the convergence threshold set to 1 × 10−5. Jackknife test was used to estimate the contribution of each environmental variable. The feature type was selected LQPTH with RM value 1. To overcome sampling bias, a bias file was created in R (v4.0.1) using the package ENMeval v0.3.0 [94]. The output format was chosen as Cloglog [95]. AUC values were examined to check for the predictive ability of the model.

Availability of data and materials

Sequences generated in this study have been deposited in GenBank under the Accession numbers OM830436–OM830486, the details of which can be found in Additional file 4: Table S1.

Change history

Abbreviations

IHR:

Indian Himalayan region

Cyt-b :

Cytochrome-b

LGM:

Last glacial maximum

IUCN:

International Union for Conservation of Nature

YBP:

Years before present

ABC:

Approximate Bayesian computation

WMP:

Western metapopulation

EMP:

Eastern metapopulation

W-S:

West Sutlej

S–B:

Sutlej_Bhagirathi

B–M:

Bhagirathi_Mahakali

M–Kr:

Mahakali_Karnali

Kr–Kg:

Karnali_KaliGandaki

Kg–My:

KaliGandaki_Marsyagandi

My-Bg:

Marsyagandi_BudhiGandaki

Bg–A:

BudhiGandaki_Arun

A–T:

Arun_Tamor

E-T:

EastTamor

HKY:

Hasegawa, Kishino and Yano

ML:

Maximum likelihood

ESS:

Effective sample size

SDM:

Species distribution modelling

AUC:

Area under curve

References

  1. Hewitt GM. Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999;68:87–112.

    Google Scholar 

  2. Geffen E, Anderson MJ, Wayne RK. Climate and habitat barriers to dispersal in the highly mobile grey wolf. Mol Ecol. 2004;13(8):2481–90.

    CAS  PubMed  Google Scholar 

  3. Thanou E, Sponza S, Nelson EJ, Perry A, Wanless S, Daunt F, et al. Genetic structure in the European endemic seabird, Phalacrocorax aristotelis, shaped by a complex interaction of historical and contemporary, physical and nonphysical drivers. Mol Ecol. 2017;26:2796–811.

    PubMed  Google Scholar 

  4. Eriksson J, Hohmann G, Boesch C, Vigilant L. Rivers influence the population genetic structure of bonobos (Pan paniscus). Mol Ecol. 2004;13(11):3425–35.

    CAS  PubMed  Google Scholar 

  5. Macfarlane CBA, Natola L, Brown MW, Burg TM. Population genetic isolation and limited connectivity in the purple finch (Haemorhous purpureus). Ecol Evol. 2016;6:8304–17.

    PubMed  PubMed Central  Google Scholar 

  6. Li R, Chen W, Tu L, Fu J. Rivers as barriers for high elevation amphibians: a phylogeographic analysis of the alpine stream frog of the Hengduan Mountains. J Zool. 2008;277:309–16.

    Google Scholar 

  7. Li Y-S, Shih K-M, Chang C-T, Chung J-D, Hwang S-Y. Testing the effect of mountain ranges as a physical barrier to current gene flow and environmentally dependent adaptive divergence in Cunninghamia konishii (Cupressaceae). Front Genet. 2019;10(742):1–15.

    PubMed  PubMed Central  Google Scholar 

  8. Sanchez-Montes G, Wang J, Arino AH, Martinez-Solano I. Mountains as barriers to gene flow in amphibians: quantifying the differential effect of a major mountain ridge on the genetic structure of four sympatric species with different life history traits. J Biogeogr. 2018;45:318–31.

    Google Scholar 

  9. Wallace AR. On the monkeys of the Amazon. Proc Zool Soc London. 1852;20:107–10.

    Google Scholar 

  10. Boubli JP, Ribas C, Lynch Alfaro JW, Alfaro ME, da Siva MNF, Pinho GM, et al. Spatial and temporal patterns of diversification on the amazon: a test of the riverine hypothesis for all diurnal primates of Rio Negro and Rio Branco in Brazil. Mol Phylogenet Evol. 2015;82(B):400–12.

    PubMed  Google Scholar 

  11. Lecompte E, Bouanani M-A, de Thoisy B, Crouau-roy B. How do rivers, geographic distance, and dispersal behavior influence genetic structure in two sympatric new world monkeys? Am J Primatol. 2017;79(7):1–12.

    Google Scholar 

  12. Naka LN, Bechtoldt CL, Henriques LMP, Brumfield RT. The role of physical barriers in the location of Avian Suture Zones in the Guiana Shield, northern Amazonia.  Am Nat. 2012;179(4):E115–32.

    PubMed  Google Scholar 

  13. Fernandes AM, Cohn-Haft M, Hrbek T, Farias IP. Rivers acting as barriers for bird dispersal in the Amazon. Rev Bras Ornitol. 2014;22(4):363–73.

    Google Scholar 

  14. Nazareno AG, Dick CW, Lohmann LG. Wide but not impermeable: testing the riverine barrier hypothesis for an Amazonian plant species. Mol Ecol. 2017;26(14):3636–48.

    CAS  PubMed  Google Scholar 

  15. Godinho MBDC, da Silva FR. The influence of riverine barriers, climate, and topography on the biogeographic regionalization of Amazonian anurans. Sci Rep. 2018;8(3427):1–11.

    Google Scholar 

  16. Oliveira EF, Martinez PA, São-pedro VA, Gehara M, Burbrink FT, Mesquita DO, et al. Climatic suitability, isolation by distance and river resistance explain genetic variation in a Brazilian whiptail lizard. Heredity (Edinb). 2018;120:251–65.

    Google Scholar 

  17. Hartmann SA, Steyer K, Kraus RHS, Segelbacher G, Nowak C. Potential barriers to gene flow in the endangered European wildcat (Felis silvestris). Conserv Genet. 2013;14:413–26.

    Google Scholar 

  18. Stojak J, Wójcik JM, Ruczy I, Searle JB, Mcdevitt AD. Contrasting and congruent patterns of genetic structuring in two Microtus vole species using museum specimens. Mammal Res. 2016;61:141–52.

    Google Scholar 

  19. Gerlach G, Musolf K. Fragmentation of landscape as a cause for genetic subdivision in bank voles. Conserv Biol. 2000;14(4):1066–74.

    Google Scholar 

  20. Mullins J, Mcdevitt AD, Kowalczyk R, Ruczynska I, Gorny M, Wojcik JM. The influence of habitat structure on genetic differentiation in red fox populations in north-eastern Poland. Acta Theriol (Warsz). 2014;59:367–76.

    Google Scholar 

  21. Milner ML, Rossetto M, Crisp MD, Weston PH. The impact of multiple biogeographic barriers and hybridization on species—level differentiation. Am J Bot. 2012;99(12):1–13.

    Google Scholar 

  22. Ayres JM, Clutton-Brock TH. River boundaries and species range size in Amazonian primates. Am Nat. 1992;140(3):531–7.

    CAS  PubMed  Google Scholar 

  23. Ribas CC, Aleixo A, Afonso CRN, Miyaki CY, Cracraft J. A palaeobiogeographic model for biotic diversification within Amazonia over the past three million years. Proc R Soc B. 2012;279:681–9.

    PubMed  Google Scholar 

  24. Kawamoto Y, Takemoto H, Higuchi S, Sakamaki T, Hart Ja, Hart TB, et al. Genetic structure of wild bonobo populations: diversity of mitochondrial DNA and geographical distribution. PLoS ONE. 2013;8(3):1–9.

    Google Scholar 

  25. Arora N, Nater A, van Schaik CP, Willems EP, van Noordwijk Ma, Goossens B, et al. Effects of Pleistocene glaciations and rivers on the population structure of Bornean orangutans (Pongo pygmaeus). Proc Natl Acad Sci USA. 2010;107(50):21376–81.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Nater A, Arora N, Greminger MP, van Schaik CP, Singleton I, Wich SA, et al. Marked population structure and recent migration in the critically endangered Sumatran orangutan (Pongo abelii). J Hered. 2013;104(1):2–13.

    PubMed  Google Scholar 

  27. Ram MS, Kittur SM, Biswas J, Nag S, Shil J, Umapathy G. Genetic diversity and structure among isolated populations of the endangered gees golden langur in Assam, India. PLoS ONE. 2016;11(8):e0161866.

    PubMed  PubMed Central  Google Scholar 

  28. Khanal L, Chalise MK, Wan T, Jiang X. Riverine barrier effects on population genetic structure of the Hanuman langur (Semnopithecus entellus) in the Nepal Himalaya. BMC Evol Biol. 2018;18(159):1–16.

    Google Scholar 

  29. Khanal L, Chalise MK, He K, Acharya BK, Kawamoto Y, Jiang X. Mitochondrial DNA analyses and ecological niche modeling reveal post-LGM expansion of the Assam macaque (Macaca assamensis) in the foothills of Nepal Himalaya. Am J Primatol. 2018;80:e22748.

    PubMed  Google Scholar 

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

    Google Scholar 

  31. Hewitt GM. Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc B Biol Sci. 2004;359(1442):183–95.

    CAS  Google Scholar 

  32. Hickerson MJ, Carstens BC, Cavender-Bares J, Crandall KA, Graham CH, Johnson JB, et al. Phylogeography’s past, present, and future: 10 years after Avise, 2000. Mol Phylogenet Evol. 2010;54(1):291–301.

    CAS  PubMed  Google Scholar 

  33. Delrieu-trottin E, Hubert N, Giles EC, Chifflet-belle P, Suwalski A, Neglia V, et al. Coping with Pleistocene climatic fluctuations: demographic responses in remote endemic reef fishes. Mol Ecol. 2020;29(12):2218–33.

    PubMed  Google Scholar 

  34. Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13.

    CAS  PubMed  Google Scholar 

  35. Tian S, Lopez-Pujol J, Wang H, Ge S, Zhang Z-Y. Molecular evidence for glacial expansion and interglacial retreat during Quaternary climatic changes in a montane temperate pine (Pinus kwangtungensis Chun ex Tsiang) in southern China. Plant Syst Evol. 2010;284:219–29.

    CAS  Google Scholar 

  36. Yan F, Zhou W, Zhao H, Yuan Z, Wang Y, Jiang K, et al. Geological events play a larger role than Pleistocene climatic fluctuations in driving the genetic structure of Quasipaa boulengeri (Anura: Dicroglossidae). Mol Ecol. 2013;22:1120–33.

    PubMed  Google Scholar 

  37. Leite YLR, Costa LP, Loss AC, Rocha RG, Batalha-filho H, Bastos AC, et al. Neotropical forest expansion during the last glacial period challenges refuge hypothesis. PNAS. 2016;113(4):1008–13.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Lessa EP, Cook JA, Patton JL. Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary. Proc Natl Acad Sci. 2003;100(18):10331–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Gutiérrez-Rodriguez J, Barbosa AM, Martinez-Solano I. Present and past climatic effects on the current distribution and genetic diversity of the Iberian spadefoot toad (Pelobates cultripes): an integrative approach. J Biogeogr. 2017;44:245–8.

    Google Scholar 

  40. Harrington SM, Hollingsworth BD, Higham TE, Reeder TW. Pleistocene climatic fluctuations drive isolation and secondary contact in the red diamond rattlesnake (Crotalus ruber) in Baja California. J Biogeogr. 2017;45(1):64–75.

    Google Scholar 

  41. Loera I, Ickert-bond SM, Sosa V. Pleistocene refugia in the Chihuahuan Desert: the phylogeographic and demographic history of the gymnosperm Ephedra compacta. J Biogeogr. 2017;44:1–11.

    Google Scholar 

  42. Horreo JL, Pelaez ML, Suarez T, Breedveld MC, Heulin B, Surget-Groba Y, et al. Phylogeography, evolutionary history and effects of glaciations in a species (Zootoca vivipara) inhabiting multiple biogeographic regions. J Biogeogr. 2018;45:1616–27.

    Google Scholar 

  43. Livraghi L, Vod R, Evans LC, Gibbs M, Dinca V, Holland PWH, et al. Historical and current patterns of gene flow in the butterfly Pararge aegeria. J Biogeogr. 2018;45:1–12.

    Google Scholar 

  44. Holder K, Montgomerie R, Friesen VL. A test of the glacial refugium hypothesis using patterns of mitochondrial and nuclear DNA sequence variation in rock ptarmigan (Lagopus mutus). Evolution. 1999;53(6):1936–50.

    CAS  PubMed  Google Scholar 

  45. Loehr J, Worley K, Grapputo A, Carey J, Veitch A, Coltman DW. Evidence for cryptic glacial refugia from North American mountain sheep mitochondrial DNA. J Evol Biol. 2006;19:419–30.

    CAS  PubMed  Google Scholar 

  46. Excoffier L, Foll M, Petit RJ. Genetic consequences of range expansions. Annu Rev Ecol Evol Syst. 2009;40:481–501.

    Google Scholar 

  47. Unmack PJ, Bagley JC, Adams M, Hammer MP, Johnson JB. Molecular phylogeny and phylogeography of the Australian freshwater fish genus Galaxiella, with an emphasis on Dwarf Galaxias (G. pusilla). PLoS ONE. 2012;7(6):e38433.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403:853–8.

    CAS  PubMed  Google Scholar 

  49. Roy AB, Purohit R. The himalayas: evolution through collision. In: Roy AB, Purohit R, editors. Indian Shield. Amsterdam: Elsevier; 2018. p. 311–27.

    Google Scholar 

  50. Srinivasan U, Tamma K, Ramakrishnan U. Past climate and species ecology drive nested species richness patterns along an east–west axis in the Himalaya. Glob Ecol Biogeogr. 2014;23(1):52–60.

    Google Scholar 

  51. Dahal N, Kumar S, Noon BR, Nayak R, Lama RP, Ramakrishnan U. The role of geography, environment, and genetic divergence on the distribution of pikas in the Himalaya. Ecol Evol. 2020;10:1539–51.

    PubMed  PubMed Central  Google Scholar 

  52. Arekar K, Sathyakumar S, Karanth KP. Integrative taxonomy confirms the species status of the Himalayan langurs, Semnopithecus schistaceus Hodgson, 1840. J Zool Syst Evol Res. 2020;59(2):543–56.

    Google Scholar 

  53. Wangchuk T, Inouye DW, Hare MP. The emergence of an endangered species: evolution and phylogeny of the Trachypithecus geei of Bhutan. Int J Primatol. 2008;29:565–82.

    Google Scholar 

  54. Groves C. Primate taxonomy. Washington, D.C: Smithsonian Institution Press; 2001.

    Google Scholar 

  55. Ross C, Srivastava A, Pirta RS. Human influences on the population density of Hanuman Langurs Presbytis entellus and Rhesus Macaques Macaca mulatta in Shimla, India. Biol Conserv. 1993;65(2):159–63.

    Google Scholar 

  56. Pirta RS, Gadgil M, Kharshikar A V. Management of the rhesus monkey Macaca mulatta and Hanuman langur Presbytis entellus in Himachal Pradesh, India. Biol Conserv. 1997;79:97–106.

    Google Scholar 

  57. Koenig A, Borries C, Chalise MK, Winkler P. Ecology, nutrition, and timing of reproductive events in an Asian primate, the Hanuman langur (Presbytis entellus). J Zool. 1997;243:215–35.

    Google Scholar 

  58. Launhardt K, Epplen C, Epplen JT, Winkler P. Amplification of microsatellites adapted from human systems in faecal DNA of wild Hanuman langurs (Presbytis entellus). Electrophoresis. 1998;19(8–9):1356–61.

    CAS  PubMed  Google Scholar 

  59. Borries C, Koenig A, Winkler P. Variation of life history traits and mating patterns in female langur monkeys (Semnopithecus entellus). Behav Ecol Sociobiol. 2001;50:391–402.

    Google Scholar 

  60. Sayers K, Norconk MA. Himalayan Semnopithecus entellus at Langtang National Park, Nepal: diet, activity patterns, and resources. Int J Primatol. 2008;29(2):509–30.

    Google Scholar 

  61. Cornuet J, Pudlo P, Veyssier J, Dehne-garcia A, Gautier M, Raphael L, et al. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics. 2014;30(8):1187–9.

    CAS  PubMed  Google Scholar 

  62. Karanth KP, Singh L, Stewart C-B. Mitochondrial and nuclear markers suggest Hanuman langur (Primates: Colobinae) polyphyly: implications for their species status. Mol Phylogenet Evol. 2010;54(2):627–33.

    CAS  PubMed  Google Scholar 

  63. Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B, et al. The last glacial maximum. Science. 2009;325:710–4.

    CAS  PubMed  Google Scholar 

  64. Roy PS, Murthy MSR, Roy A, Kushwaha SPS, Singh S, Jha CS, et al. Forest fragmentation in India. Curr Sci. 2013;105(6):774–80.

    Google Scholar 

  65. Bagaria P, Sharma LK, Joshi BD, Kumar H, Mukherjee T, Thakur M, et al. West to east shift in range predicted for Himalayan Langur in climate change scenario. Glob Ecol Conserv. 2020;22:e00926.

    Google Scholar 

  66. Yang M, Yang Y, Cui D, Fickenscher G, Zinner D, Roos C, et al. Population genetic structure of Guizhou snub-nosed monkeys (Rhinopithecus brelichi) as inferred from mitochondrial control region sequences, and comparison with R. roxellana and R. bieti. Am J Phys Anthropol. 2012;147:1–10.

    CAS  PubMed  Google Scholar 

  67. Slatkin M, Hudson RR. Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991;129:555–62.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Harcourt AH, Wood MA. Rivers as barriers to primate distributions in Africa. Int J Primatol. 2012;33(1):168–83.

    Google Scholar 

  69. Hugall A, Moritz C, Moussalli A, Stanisic J. Reconciling paleodistribution models and comparative phylogeography in the wet tropics rainforest land snail Gnarosophia bellendenkerensis (Brazier 1875). Proc Natl Acad Sci. 2002;99(9):6112–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. Gaubert P, Njiokou F, Ngua G, Afiademanyo K, Dufour S, Malekani J, et al. Phylogeography of the heavily poached African common pangolin (Pholidota, Manis tricuspis) reveals six cryptic lineages as traceable signatures of Pleistocene diversification. Mol Ecol. 2016;25:5975–93.

    CAS  PubMed  Google Scholar 

  71. Schmidt D, Pool J. The effect of population history on the distribution of the Tajima’s D statistic. New York: Cornell University Press; 2002.

    Google Scholar 

  72. Liu Z, Ren B, Wei F, Long Y, Hao Y, Li M. Phylogeography and population structure of the Yunnan snub-nosed monkey (Rhinopithecus bieti) inferred from mitochondrial control region DNA sequence analysis. Mol Ecol. 2007;16(16):3334–49.

    CAS  PubMed  Google Scholar 

  73. Li M, Liu Z, Gou J, Ren B, Pan R, Su Y, et al. Phylogeography and population structure of the golden monkeys (Rhinopithecus roxellana): inferred from mitochondrial DNA sequences. Am J Primatol. 2007;69:1195–209.

    CAS  PubMed  Google Scholar 

  74. Benn DI, Owen LA. The role of the Indian summer monsoon and the mid-latitude westerlies in Himalayan glaciation: review and speculative discussion. J Geol Soc. 1998;155(2):353–63.

    Google Scholar 

  75. Owen LA, Finkel RC, Caffee MW. A note on the extent of glaciation throughout the Himalaya during the global last glacial maximum. Quat Sci Rev. 2002;21(1–3):147–57.

    Google Scholar 

  76. Eugster P, Scherler D, Thiede RC, Codilean AT, Strecker MR. Rapid last glacial maximum deglaciation in the Indian Himalaya coeval with midlatitude glaciers: new insights from 10Be-dating of ice-polished bedrock surfaces in the Chandra Valley, NW Himalaya. Geophys Res Lett. 2016;43(4):1589–97.

    Google Scholar 

  77. Edgar RC. Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32(5):1792–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  78. Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33(7):1870–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  79. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.

    CAS  PubMed  PubMed Central  Google Scholar 

  80. Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22(21):2688–90.

    CAS  PubMed  Google Scholar 

  81. Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A, Höhna S, et al. Mrbayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61(3):539–42.

    PubMed  PubMed Central  Google Scholar 

  82. Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. 2013. http://beast.bio.ed.ac.uk/Tracer.

  83. Bouckaert R, Vaughan TG, Barido-Sottani J, Duchene S, Fourment M, Gavryushkina A, et al. BEAST 2.5: an advanced software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2019;15(4):e1006650.

    CAS  PubMed  PubMed Central  Google Scholar 

  84. Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29(6):1695–701.

    CAS  PubMed  Google Scholar 

  85. Soares P, Ermini L, Thomson N, Mormina M, Rito T, Rohl A, et al. Correcting for purifying selection: an improved human mitochondrial molecular clock. Am J Hum Genet. 2009;84:740–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  86. Rozas J, Ferrer-mata A, Sanchez-DelBarrio JC, Guirao-rico S, Librado P, Ramos-onsins SE, et al. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34(12):3299–302.

    CAS  PubMed  Google Scholar 

  87. Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16(1):37–48.

    CAS  PubMed  Google Scholar 

  88. Leigh JW, Bryant D. POPART: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6:1110–6.

    Google Scholar 

  89. Barido-Sottani J, Boskova V, du Plessis L, Kuhnert D, Magnus C, Mitov V, et al. Taming the BEAST—a community teaching material resource for BEAST 2. Syst Biol. 2018;67(1):170–4.

    PubMed  Google Scholar 

  90. Fagundes NJR, Ray N, Beaumont M, Neuenschwander S, Salzano FM, Bonatto SL, et al. Statistical evaluation of alternative models of human evolution. PNAS. 2007;104(45):17614–17619.

    CAS  PubMed  PubMed Central  Google Scholar 

  91. Minhas RA, Ahmed KB, Awan MS, Zaman Q, Dar NI, Ali H. Distribution patterns and population status of the Himalayan grey langur (Semnopithecus ajax) in Machiara National Park, Azad Jammu and Kashmir, Pakistan. Pak J Zool. 2012;44(3):869–77.

    Google Scholar 

  92. Minhas RA, Khan MN, Awan MS, Ahmad B, Bibi SS, Hanif M, et al. RAPD based genetic diversity of endangered Himalayan Gray Langur (Semnopithecus ajax) populations of Pakistan. Pak J Zool. 2018;50(6):2059–71.

    CAS  Google Scholar 

  93. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell. 2006;190:231–59.

    Google Scholar 

  94. Muscarella R, Galante PJ, Soley-guardia M, Boria RA, Kass JM, Uriarte M, et al. ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for MAXENT ecological niche models. Methods Ecol Evol. 2014;5:1198–205.

    Google Scholar 

  95. Phillips SJ, Anderson RP, Dudík M, Schapire RE, Blair ME. Opening the black box: an open-source release of Maxent. Ecography. 2017;40:887–93.

    Google Scholar 

Download references

Acknowledgements

We would like to thank the DBT-IISc partnership program for providing funds for laboratory work. We would like to thank Dr. Chetan Nag and Dr. Maitreya Sil for providing some of the samples used in this study. Dr. Bilal Habib for help during sample collection We would also like to thank the forest departments of Uttarakhand, Himachal Pradesh and Sikkim for providing necessary collection permits. Additionally, thanks to the Department of Wildlife Protection, Jammu and Kashmir for permits provided to MK. KA would like to thank the Rufford foundation and Primate Conservation Inc. for providing funding for field work. KA would also like to thank ANCF for the administrative work in managing the finances of the project. MK would like to thank Tahir Gazanfar for samples collected from Limber. KA would like to thank Sartaj, Suresh, Rajat, Bhim Singh, Kunzang and the Bhutia family for their help during fieldwork in different parts of the Himalayas.  We also thank the two anonymous reviewers for their comments which significantly improved the manuscript.

Funding

This work was supported by the Ruffords small grant for field work (18855-1) and a fieldwork grant from Primate Conservation Inc. (001452) to Kunal Arekar. The molecular lab work and data generation was undertaken with support from the annual DBT-IISc partnership fund to Praveen Karanth.

Author information

Authors and Affiliations

  1. Centre for Ecological Sciences, Indian Institute of Science, Bengaluru, India

    Kunal Arekar, Neha Tiwari & Praveen Karanth

  2. Wildlife Institute of India, Dehradun, India

    Kunal Arekar & Sambandam Sathyakumar

  3. Wildlife Research and Conservation Foundation, Srinagar, Jammu and Kashmir, India

    Mehreen Khaleel

Authors
  1. Kunal Arekar
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Neha Tiwari
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Sambandam Sathyakumar
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Mehreen Khaleel
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Praveen Karanth
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

KA, PK and SS conceptualised the study. KA collected the samples, generated some of the molecular sequences, analysed all the data, and prepared all the figures and tables. NT generated some of the molecular sequences. MK collected the samples from Jammu and Kashmir. KA and PK wrote the paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Kunal Arekar.

Ethics declarations

Ethics approval and consent to participate

The ethics approval was not needed since we collected faecal samples. The field work was carried out in different states under the following permits from the respective forest departments—Jammu and Kashmir (WLP/Res/16-17/338-44), Himachal Pradesh (WL/Research-Study/2211), Uttarakhand (599/5-6), Sikkim (67/GOS/FEWMD/BD-R-2014/CCF (T&HQ) 2/Pro/Res).

Consent of publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

The original online version of this article was revised: Following the publication of the original article, we were notified that the order of the first and last name of the 3rd author was reversed. “Sathyakumar Sambandam” should be “Sambandam Sathyakumar”.

Supplementary Information

Additional file 1: Figure S1.

Maximum likelihood(ML) phylogeny of Himalayan langur for mitochondrial cytochrome b (Cyt-b) gene.Numbers at the node indicate ML bootstrap (MLBS) values. Node support values of > 75 are shown. Node support values only for major clades are shown.

Additional file 2: Figure S2.

Divergence time tree for the mitochondrial cytochrome b (Cyt-b) dataset constructed in BEASTv2.6. The bars at the nodes indicate the 95% credible interval. Details of the collapsed IHR node can be found in Additional file 4: Table S3. IHR = Indian HimalayanRegion; WMP = Western metapopulation; EMP = Eastern metapopulation.

Additional file 3:Figure S3.

Response curves for top three variables of importance for Himalayan langur in the niche modelling analysis.

Additional file 4: Table S1.

Samples used for molecular phylogenetic analysis. Samples that are marked with † were sequenced in this study. Samples with the same accession numbers are identical sequences. Table S2. Sequences used for divergence dating analysis in BEAST v2.6.6. Sequences that are marked with † were sequenced in this study. Sequences from sr. no. 1–66 constitutes the collapsed clade IHR in Additional file 2: Fig. S2. Table S3. Model selection for Maxent analysis—the table shows AUC values for different models. AUC values in bold shows the features and RM values selected.

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

Verify currency and authenticity via CrossMark

Cite this article

Arekar, K., Tiwari, N., Sathyakumar, S. et al. Geography vs. past climate: the drivers of population genetic structure of the Himalayan langur. BMC Ecol Evo 22, 100 (2022). https://doi.org/10.1186/s12862-022-02054-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-022-02054-1

Keywords

  • Genetic diversity
  • Phylogeography
  • Riverine barrier effect
  • Past climate
  • Heterogenous landscape

BMC Ecology and Evolution

ISSN: 2730-7182

Contact us