royalsocietypublishing.org/doi/full/10.1098/rspb.2017.1804 Evolutionary history of enigmatic bears in the Tibetan Plateau–Himalaya region and the identity of the yeti.
Although anecdotally associated with local bears (Ursus arctos and U. thibetanus), the exact identity of ‘hominid’-like creatures important to folklore and mythology in the Tibetan Plateau–Himalaya region is still surrounded by mystery. Recently, two purported yeti samples from the Himalayas showed genetic affinity with an ancient polar bear, suggesting they may be from previously unrecognized, possibly hybrid, bear species, but this preliminary finding has been under question. We conducted a comprehensive genetic survey of field-collected and museum specimens to explore their identity and ultimately infer the evolutionary history of bears in the region. Phylogenetic analyses of mitochondrial DNA sequences determined clade affinities of the purported yeti samples in this study, strongly supporting the biological basis of the yeti legend to be local, extant bears. Complete mitochondrial genomes were assembled for Himalayan brown bear (U. a. isabellinus) and black bear (U. t. laniger) for the first time. Our results demonstrate that the Himalayan brown bear is one of the first-branching clades within the brown bear lineage, while Tibetan brown bears diverged much later. The estimated times of divergence of the Tibetan Plateau and Himalayan bear lineages overlap with Middle to Late Pleistocene glaciation events, suggesting that extant bears in the region are likely descendants of populations that survived in local refugia during the Pleistocene glaciations.
The Tibetan Plateau, the most extensive and highest plateau in the world with an average altitude of 4500 m above sea level, is partly surrounded by the Himalayan range and many of Earth's highest mountains. Dramatic environmental changes caused by the uplift of the plateau and climatic oscillations during the Quaternary glaciations substantially impacted the evolution, diversification, and distribution of local plant and animal species [1]. Because of its heterogeneous habitat and topography, the region sustains a distinct biome with rich biological diversity and high level of endemism [2]. Extant plants and animals on the plateau are likely either descendants of relict colonists that migrated from other areas or recently derived endemic species [3–10]. However, the colonization and population expansion history of many species remains poorly understood, despite current and future impacts of climate change and anthropogenic threats to diversity loss.
Two brown bear subspecies, the Himalayan (Ursus arctos isabellinus) and the Tibetan (U. a. pruinosus) brown bear, inhabit the northwestern Himalayan region and southeastern Tibetan Plateau, respectively [11–14] (figure 1). These two subspecies have distinct skull features and the Himalayan brown bear is characterized by its paler and reddish-brown fur, while the Tibetan brown bear has generally darker fur with a developed, white ‘collar’ around the neck [11]. As the most widely distributed bear in the world, phylogeography of the brown bear has been well studied in North America, Europe and Japan [10,16–24]. However, due to limited sampling, very few studies have been conducted on these enigmatic subspecies. Two complete mitochondrial genomes (mitogenomes) from captive Tibetan brown bears are available, while only two short fragments of mitochondrial DNA (mtDNA) sequences from the Himalayan brown bear have been published [10,15]. Phylogenetic analyses based on these sequences suggested that the Tibetan brown bear might be a relict population of the Eurasian brown bear [10], and that the Himalayan brown bear, which is genetically distinct from the Tibetan brown bear, may represent a more ancient lineage [15]. However, phylogenetic relationships deduced from limited genetic data and number of individuals have put these preliminary findings into question. For example, the phylogenetic placement of a Gobi brown bear (U. a. gobiensis) sequence [25] was inconsistent with a later study also including sequences from Himalayan brown bear [15], and phylogenetic trees based on mtDNA control region and cytochrome b sequences, respectively, of the Tibetan brown bear were incongruent [26]. The other bear species found to inhabit the Tibetan Plateau–Himalaya region is the Asian black bear (U. thibetanus), which historically had a continuous distribution from southeastern Iran through Afghanistan and Pakistan to India, Nepal, China, Korea, Japan, and south into Myanmar and the Malayan peninsula [12,27,28]. Today it occupies a patchy distribution throughout its historic range, including across a narrow band from Pakistan, Kashmir and to Bhutan, the home range of the Himalayan black bear (U. t. laniger) [27,29], which was described as distinguished from other black bear populations by its longer, thicker fur and smaller, whiter chest mark [11]. Although the range of Asian black bear overlaps with brown bear in the Tibetan Plateau–Himalaya region, it is mostly found at lower altitudes in forested hills ranging from 1200 to 3300 m [12,29]. So far, little is known about the evolutionary history of black bear in the region and no sequence data are available from the Himalayan black bear. To elucidate the evolutionary and migration history of the Himalayan and Tibetan bears, more genetic data from additional individuals are critically needed.
It has been reported that the brown bear populations in the Tibetan Plateau–Himalaya region have declined by more than half in the past century because of habitat loss, fragmentation, poaching and intense hunting by humans [12,29–31]. Facing the same threats as brown bears, Asian black bear populations have also decreased in the past few decades [29,32,33]. The Himalayan brown bear is listed in the IUCN (International Union for the Conservation of Nature) red list of threatened species as critically endangered [34], while the Asian black bear is listed as vulnerable [27]. Hence, clarifying population structure and genetic diversity for conservation management purposes is also urgently needed for these endangered bear species.
The Tibetan Plateau–Himalaya region is also known for the legend of purported ‘hominid’-like creatures, referred to as the ‘yeti’, ‘chemo’, ‘mheti’ or ‘bharmando’, among other regional monikers (for simplicity they are referred to in this paper as yeti). Despite decades of research and anecdotal association with bears and other mammals in the region [35,36], the species identity of the mysterious yeti is still debated, given the lack of conclusive evidence. A survey of hair samples attributed to yeti and other anomalous, supposed primates, was recently conducted to identify their genetic affinities [37]. Based on a short fragment of the mtDNA 12S rRNA gene from two samples collected in Ladakh, India and Bhutan, respectively, and a 100% match to a sequence recovered from a subfossil polar bear [38], Sykes et al. [37] speculated that an unclassified bear species or hybrid of polar bear and brown bear might be present in the Tibetan Plateau–Himalaya region. However, this speculation was critiqued by others [39,40], and their phylogenetic analyses using the sequences from Sykes et al. and other available Ursidae sequences did not rule out the possibility that the samples belonged to brown bear. Thus, to get accurate species identification, comprehensive phylogenetic analyses using genetic information from more variable and informative loci are needed.
Here, we report on new analyses of 24 field-collected and museum specimens, including hair, bone, skin and faecal samples, collected from bears or purported yetis in the Tibetan Plateau–Himalaya region. Based on both amplified mtDNA loci as well as complete mitogenomes, we reconstructed maternal phylogenies to increase knowledge about the phylogenetic relationships and evolutionary history of Himalayan and Tibetan bears.
A total of 24 samples, including hair, tissue, bone and faeces, were analysed in this study (electronic supplementary material, table S1). Of these, 12 samples had been collected for a previous analysis of Himalayan brown bear in the Khunjerab National Park, Northern Pakistan [30], two samples were from purported Himalayan brown bears housed in the Lahore and Islamabad Zoos, one bone sample (M-70448) recorded as U. a. pruinosus was obtained from the American Museum of Natural History, and nine samples were provided to us by the Reinhold Messner Museum and the Icon Film Company.
Genomic DNA from 12 faecal samples collected in the Khunjerab National Park, Northern Pakistan [41], were previously extracted using the QIAmp DNA Stool Kit (Qiagen, USA) in a room dedicated to processing hairs and faeces [30]. DNA from two ethanol-preserved hair samples from Lahore and Islamabad Zoos were isolated in a room dedicated to nucleic acid extraction from modern samples. A DNeasy Blood & Tissue DNA Kit (Qiagen, USA) was used according to the manufacturer's protocol, except for the following modifications to optimize extraction of DNA from hair: 10 strands of hair from each sample were cut into fragments of approximately 0.5 cm with a sterile razor blade. Ethanol was allowed to evaporate (approx. 1 h), and hair fragments were transferred to a microcentrifuge tube. Three hundred microlitres of ATL buffer, 20 µl proteinase K, 20 µl 1M DTT (dithiothreitol) and 4 µl RNase A were added, and samples were incubated at 56°C overnight until completely lysed. A negative control was prepared alongside each hair sample. Following lysis, 300 µl AL buffer and 300 µl 100% EtOH were added to each sample, and the mixture was pipetted into the DNeasy Mini Spin Column and centrifuged for 2 min. DNA was eluted twice with 50 µl AE buffer for a total elution volume of 100 µl. The remaining 10 samples, which had not been intentionally preserved for later extraction of DNA, were regarded as non-modern (ancient) samples, and thus DNA extractions and pre-amplifications were performed in a dedicated state-of-the-art cleanroom facility, physically separated from any modern DNA laboratory and appropriate for ancient DNA research. The following protocols designed for ancient DNA extraction were used: for bone samples, 50–100 mg fine bone powder was obtained from each sample by using a dental drill (HKM Surgical Handpiece, Pearson Dental, USA), and 50–100 mg skin samples were sliced into approximately 1 mm pieces with a sterile razor blade. DNA from the bone powder and the sliced skin samples was extracted using the protocol in Dabney et al. [42]. DNA from the hair samples were extracted using the protocol provided by Gilbert et al. [43] with the following modifications: 1 ml digestion buffer was used for each hair extraction. After purification with phenol and chloroform, additional purification was performed using Qiagen MinElute PCR Purification Kit (Qiagen, USA). Finally, a 12.5 µl EB buffer elution step was performed twice to obtain a total elution volume of 25 µl. DNA from approximately 100 mg faecal samples was extracted using the QIAmp DNA Stool Kit (Qiagen, USA). The final elution step was also performed twice to obtain a total volume of 100 µl. Negative controls were prepared alongside all extractions.
PCR amplifications from modern DNA were performed in a 25 µl reaction volume each containing 2.5 µl of 10 × PCR buffer (Applied Biosystems, USA), 1.0 µl of dNTP mixture (2.5 mM each dNTP; Applied Biosystems), 2.5 µl of MgCl2 (25 mM, Applied Biosystems), 0.1 µl of Taq DNA polymerase (5–10 U µl−1; Applied Biosystems, AmpliTaq Gold), 1 µl each of the forward and reverse primers (10 µM), 2 µl of the genomic DNA and 17.4 µl of H2O. The PCR reaction mix for ancient DNAs was prepared in the cleanroom by adding 21 µl H2O, 1 µl of each forward and reverse primer, and 2 µl genomic DNA to each GE illustra PuReTaq Ready-To-Go PCR bead (GE Healthcare, USA). A touchdown thermal cycling protocol was used as follows: 10 min at 94°C, 10 cycles of 30 s at 94°C, 30 s annealing with the temperature decreasing every cycle by 0.5°C from 55°C to 50°C, and 30 s extension at 72°C, followed by 25 cycles the annealing temperature set to 50°C and denaturation and extension phases as above. For samples of unknown identity, two sets of mtDNA 12S rRNA primers [44,45] were used to determine their approximate taxonomic affinity. Bear-specific primers targeting the mtDNA control region and cytochrome b ([46] and primers designed for this study; see electronic supplementary material, table S2) were used for samples identified as ursid bears. PCR products were Sanger sequenced directly using the same primers as in the PCR.
Fifty microlitres of DNA extracts from four samples were sent to MYcroarray (http://www.mycroarray.com) for preparation of Ion Torrent sequencing libraries and mtDNA target enrichment and sequencing, using the following protocol. Sample libraries were quantified using spectrofluorometry, which indicated between 5 and 255 total nanograms (0.2–8.5 ng µl−1) of double-stranded DNA. Each library was then individually target enriched using a custom-designed ursid mitogenome bait set manufactured by MYcroarray. The standard MYbaits v. 3.0 protocol was applied with hybridization for 21 h at 60°C at all relevant steps. Following clean up, half of each bead-bound library was amplified in a 50 µl reaction with universal Ion Torrent adapter-primers for 10 cycles using a KAPA HiFi premix (KAPA Biosystems) and the manufacturer's recommended thermal profile coupled with 62°C annealing temperature. After amplification, the beads were pelleted and the supernatant was purified using SPRI beads and eluted in Tris-HCl buffer containing 0.05% Tween-20. The enriched libraries were quantified with spectrofluorometry, which indicated between 1.12 and 4.21 total nanograms dsDNA per library (0.03–0.12 ng µl−1). Equal masses of each library were pooled, bead-templated and sequenced alongside other project libraries on the Ion Proton platform using the Ion PI Chip Kit v2 chemistry. Following sequencing, reads were de-multiplexed, quality trimmed and filtered using the default settings on the Ion Torrent Suite v. 4.4.3.
Assembly of mitochondrial genomes was performed using the following strategy: species-specific mitochondrial reference genomes were selected from initial species identification based on phylogenetic analyses of amplicon sequences (results not shown). All Ion Torrent reads were first aligned against the reference genome using BWA aln (v. 0.7.13) [47] using the default parameters, except for the parameter ‘-l 1024’ to disable the seed and increase high-quality hits for the damaged ancient DNA reads [48]. The remaining unmapped reads were then aligned against the same reference using BWA mem with default parameters (see electronic supplementary material, table S3, for assembly statistics). We filtered for human contamination by applying an edit-distance based strategy [48]. All reads were mapped to a human mitochondrial genome reference (NCBI accession J01415.2) using the same BWA mapping method described above. Reads with a higher mapping edit-distance to human mtDNA than to bear mitochondrial genomes were considered of likely human origin and were removed from the bear mitogenome mapping results. PCR duplicates were removed with the MarkDuplicates tool in the Picard software suite v. 1.112 using lenient validation stringency (http://broadinstitute.github.io/picard/). Consensus calling was carried out using Samtools mpileup [49] with default settings.
Complete mitochondrial genomes, partial control region sequences, and cytochrome b sequences for 11 Asian black bears, 76 American black bears, two cave bears (U. spelaeus), 200 brown bears, and 52 polar bears were obtained from GenBank (electronic supplementary material, table S4). Two GenBank datasets were created: one dataset included only complete mitogenomes for the non-Tibetan/Himalayan bears and both partial (amplicon sequences) and complete mitogenomes for Tibetan and Himalayan bear lineages, while the other dataset included both amplicon sequences and complete mitogenomes for non-Tibetan/Himalayan bears. All new sequences produced in this study were added to these two GenBank datasets and used in the phylogenetic analyses. Sloth bear (U. ursinus) and sun bear (U. malayanus) sequences were included to root the trees (electronic supplementary material, table S4). Alignments were generated using MAFFT [50] followed by manual adjustment in BioEdit [51] to exclude the variable number tandem repeats of the D-loop. The total length of the final alignment was 16 412 bp. Maximum-likelihood (ML) phylogenetic analyses were performed using RAxML-HPC BlackBox v. 8.2.8 [52] in the CIPRES Science Gateway under the GTR substitution model, which was identified as the best-supported model by jmodeltest2 [53,54]. A total of 1000 bootstrap replicates were conducted to evaluate branch support. Bayesian inference (BI) phylogenetic analyses were carried out using MrBayes v. 3.2.6 [55] in two runs of 5 000 000 Markov chain Monte Carlo (MCMC) generations, with trees for estimation of the posterior probability distribution sampled every 100 generations. The best-fit substitution model was determined by the program by setting Nst=mixed; 500 000 trees were discarded as burn-in.
Much more info on site; but I feel that this is a huge overload already. So why this post? Two very different brown bear subspecies live at higher altitudes than any other bear species - with the Andes bear running a close second place.
