Phylogenetic analysis of marginal Asiatic black bears reveals a recent Iranian–Himalayan divergence and has implications for taxonomy and conservation

A small population of Asiatic black bear—known as the Baluchistan black bear—survives in the western limit of the species’ range in Iran, where the species is rare, difficult to monitor and occupy an atypical habitat with extreme environmental conditions. Through the use of noninvasively collected samples, we analyzed mitochondrial DNA control region sequences to evaluate the phylogenetic relationships and divergence time between the Baluchistan Iranian black bear population and other Asian populations. Phylogenetic analyses indicate that Baluchistan and Nepalese (Himalayan) populations are monophyletic, with their divergence time estimated at circa 120 thousand years ago. The results reveal the low level of mitochondrial DNA variability in this small and marginal population, as is the case for many bear populations living in areas with similar conditions. The divergence time between the populations from Iran and Nepal dates to the Late Pleistocene, pointing to a transitional period between colder (glacial) and warmer (interglacial) conditions that allowed forests to expand and opened new habitats to population expansions. Pending further genetic and morphological corroboration, these preliminary results suggest that all Baluchistan and Himalayan (Nepalese) black bears might be considered as synonymous under the prior U. t. thibetanus trinomial (with gedrosianus just as junior synonym). Conservation efforts on this small and endangered population remain poor, and further measures are required to guarantee its long-term survival in Iran.


Introduction
The Asiatic black bear Ursus thibetanus Cuvier (1823) is widely distributed throughout Asia, ranging from the Manchuria Region in eastern Russia to Iran, in the southwestern part of this vast continent (Garshelis and Steinmetz 2016). Unlike their eastern Asian counterparts, the populations whose distribution area spans from China to Iran remain poorly studied. As such, it is not surprising that the genetic status of these populations is still hardly known; in contrast, genetic tools have been widely employed to study the eastern Asian populations (e.g., Saitoh et al. 2001;Ishibashi and Saitoh 2004;Ohnishi et al. 2007Ohnishi et al. , 2009Yasukochi et al. 2009;Choi et al. 2010;Kim et al. 2011;Uno et al. 2015;Wu et al. 2015). To date, only one study (Kadariya et al. 2018) has employed molecular markers to assess Asiatic black bear genetic diversity east from China. Interestingly, this study showed that the Nepalese Asiatic black bear population was in fact a distinct lineage and basal to all other mainland populations.
In Iran, Asiatic black bears are confined to small and isolated patches in the remote mountainous areas of the southeastern part of the country (Almasieh et al. 2016;Almasieh and Kaboli 2018;Farashi and Erfani 2018). This population, also known as the Baluchistan black bear, is the westernmost marginal population of the species, as well as one of the most threatened (Fahimi et al. 2011;Garshelis and Steinmetz 2016). Furthermore, not only is this population geographically marginal but it is also an ecologically marginal population. In the dry landscape and sparse woodlands of this area that are intersected by lowland deserts and wide plateaus, the bear inhabits a suboptimal habitat with extreme environment conditions that differs from the other parts of species' range in Asia (Yusefi 2013). Currently, this population is under severe threat, primarily by habitat loss and fragmentation, along with human persecution by human-bear conflicts (Ghadirian et al. 2017;Fahimi et al. 2018) and, therefore, it is listed as "Endangered" at a national level (Yusefi et al. 2019). Due to the scarcity of suitable habitat, it is believed that the Baluchistan black bear population may occur at very low densities with a population under 200 individuals, divided in scattered subpopulations, each only harboring a few individuals (Yusefi 2013). Despite this, the Baluchistan black bear is one of the least-known large mammals in the country and information and scientific data remains deficient, with a subsequent lack of efficient conservation measures (Yusefi 2013;Almasieh and Kaboli 2018).
The Baluchistan black bear was long considered as a distinct subspecies, U. t. gedrosianus (Blanford 1877). However, the taxonomic validity of U. t. gedrosianus is uncertain and the affinity of this population with that of other populations has not been tested with molecular methods (Kadariya et al. 2018). The latter can be particularly important since the taxonomy of Asiatic black bears below the species level remains uncertain, and there is no consensus with regard to either their systematics or their taxonomic nomenclature (Choi et al. 2010;Wu et al. 2015). In this study, we use noninvasively collected samples from Iran and published data from GenBank, to explore the phylogenetic relationship and divergence time of Baluchistan black bears in relation to other Asiatic black bear subspecies and provide insights into the subspecific taxonomy of the species throughout its range in Asia. To foresee management decisions and conservation planning there is a pressing need to assess the genetic status of small, fragmented and threatened populations.

Study sites, sampling and DNA extraction
Baluchistan black bears were noninvasively sampled by collecting scats in six geographically separated areas: (1) Bahr-e Asman, (2) Dehbakri-Dalfard (Zaryab), (3) Kahnouj County and (4)   (2), Kahnouj County (3), Marz-Bashagard (4), Roudan County (5), and Nikshahr County (6) were sampled, with the exception of Jebal Barez (7) 2013). Although the probable range of the species is wide, its presence has merely been confirmed in some localities in southeast Iran, and all these areas with the exception of Jebal Barez were sampled (Fig. 1). Overall, the average distance between localities was circa 180 km (distance; minimum 40 km, maximum 500 km). 218 fecal samples were collected opportunistically between 2008 and 2011 along animal trails, or from cave entrances and orchards, but only those believed to be fresh (n = 74) were used for genetic analysis. Besides the feces samples, two tissue samples belonging to dead individuals (from Dehbakri-Dalfard and Roudan) were also included in this study. Finally, one fresh scat sample from a captive bear (originally from Nikshahr area) was also included. Air-dried samples were stored inside a Ziploc plastic bag containing silica gel or in a container with pure ethanol until transferring to the Research Center in Biodiversity and Genetic Resources (CIBIO) (Porto, Portugal).

DNA extraction and sequencing
DNA extraction followed Costa et al. (2017) with overnight digestion with proteinase K and carried on a laminar flux chamber, physically isolated from the PCR room. The samples were processed in batches with a maximum of 16 samples per set. All the materials used for the extraction process was disinfected between samples and between all groups of samples as well. Negative controls were used to assess for contamination with each extraction batch.
We sequenced a 338 bp fragment of the mtDNA control region (CR) with primers UTCR1F (5′-CCT AAG ACT AAG GAA GAA G-3′) and UTCR2R (5′-TAC TCG CAA GGA TTG CTG G-3′) (Miller et al. 2006). This short sequence was targeted due to bad preservation of samples and because it holds the highest variability of the CR (Yusefi 2013). In addition, two individuals (see results) were amplified for 680 bp (which included the shorter 338 bp) with the same forward primer and the reverse primer (5′-CGT GTT CYY CGA TTC AGT GGT ATT -3′, this study) to include in the phylogenetic analyses. PCR conditions were as described in Yusefi (2013). Thirty three sequences, which represent 19 individuals identified by probability of identity through microsatellite genotyping (Yusefi 2013) were used in this study.

Data analyses
Sequences were inspected and edited using Sequencer v5.4.6 (Gene Codes Corporation, Ann Arbor, MI, USA) and thereafter submitted to GenBank (under accession numbers MT432185-MT432186). Subspecies-level identification was performed using GenBank's BLASTn search (Altschul et al. 1990). Seaview v.4.2.11 (Gouy et al. 2010) was used for preliminary alignments of sequences and were aligned thereafter in MAFFT (Katoh et al. 2002), and phylogenetic analyses were conducted using 680 bp alignment of the control region. The most appropriate substitution model (KKY + I) was determined according to the Akaike's Information Criterion (Akaike 1974) in jModeltest 2 (Posada 2008).

Phylogenetic analysis and divergence time
Phylogenetic analyses were conducted using Bayesian Inference (BI) in MrBayes v3.2.6 (Ronquist et al. 2012) and Maximum Likelihood (ML) in Garli v2.0.1 (Zwickl 2006). MrBayes was used with default priors and Markov chain settings, and with random starting trees. Each run consisted of four chains of 10 million generations, sampled every 10,000 generations. ML analysis was performed using 10 independent searches and 1000 bootstrap replicates. PAUP (Swofford 2002) was used to summarize non-parametric bootstrap support values for the best tree, after generating a majority-rule consensus tree. All analyses were performed using the CIP-RES platform (Miller et al. 2010). Consensus tree inferred for each phylogenetic approach was visualized and rooted using FigTree v1.4.3 (Rambaut 2017).
We performed statistical coalescent and phylogeographic analyses using BEAST v1.8.3 (Drummond et al. 2012; https ://beast .bio.ed.ac.uk). We inferred a mitochondrial gene tree based on a constant coalescent model under the HKY + I model strict clock. To time-calibrate the population tree, we fixed the mutation rate in the control region to 2.9 × 10 -8 substitutions/site/year following Wu et al. (2015). This data set was much reduced to include only representative of all Asiatic black bears subspecies from Japan, Korea, Russia, Nepal, Taiwan, and China (downloaded from GenBank; Supplementary material, Table S1) and combined with haplotypes from Iran (this study). We ran two independent MCMC chains, each with 10 million generations and sampling every 10,000 generations. Independent runs were evaluated for convergence and mixing by observing and comparing traces of each statistic and parameter in Tracer v1.7 (Rambaut et al. 2018; https ://beast .bio.ed.ac.uk/trace r). We considered effective sampling size (ESS) values > 200 to be good indicators of parameter mixing. The first 20% generations of each run were discarded as burnin, and samples were merged using LogCombiner v1.8.3. The resulting chains were summarized using TreeAnnotator v1.8.3, where a maximum-clade-credibility (MCC) tree with mean values were generated under the "-heights ca" option (Heled and Bouckaert 2013).
Two median-joining networks were constructed in Pop-ART, Bandelt et al. 1999; https ://popar t.otago .ac.nz) to assess the phylogenetic relationships of all known control region sequences available for Asiatic black bears (including translocated animals from captive centers and zoo animals 1 3 with uncertain origin) and a smaller network with only animals from known origin (downloaded from GenBank; Supplementary material, Table S1). We include only two U. t. japonicus haplotypes to show the position of this highly divergent clade relative to the other subspecies. The larger and smaller network alignments had 121 and 88 sequences, respectively. However, total exclusion of sequences with considerable missing data resulted in 85 and 54 sequences, respectively.

Results
Two haplotypes differing by an indel were recovered from all Iranian samples (Supplementary material, Table S2). The longer control region sequences (680 bp) belonged to each of the two haplotypes, and despite a much increase sequenced length, they differed only by the same indel. The control region belonging to a Chinese specimen (GenBank accession No. DQ402478, the complete mitochondrial genome) was highly divergent to all mainland U. thibetanus subspecies and even more divergent than U. t. japonicus. The GenBank blast searchers of the control region sequence matched U. americanus (Pallas 1780) with 98% identity (as it has been reported previously by Wu et al. 2015), despite the authors reporting a 90% similarity to the complete control region (1422 bp) of the American black bear (Hou et al. 2007). The origin on this sequence derives from Asiatic black bears collected from Sichuan Province Traditional Chinese Medicinal Materials Company and, therefore, may suggest either an unknown sample or cross contamination. Following these results, we did not name this subspecies in our analyses U. t. mupinensis (unlike analyses including its complete mitochondrial genome, Wu et al. 2015;Kadariya et al. 2018) and was excluded from the phylogenetic analyses.
The large network recovered Taiwanese, Chinese, Vietnamese, and Korean samples scattered throughout, with no evident structure or distribution patterns in relation to current distribution suggesting mixture of Asiatic black bears in recovery centers and zoos (Supplementary material, Fig. S1). The smaller network that excluded all those likely mixed individuals with uncertainty in the locality of origin (from zoos, recovery centers, or introduced animals) recovered five well-distinguishable U. thibetanus clusters (subspecies groups), in accordance with Kadariya et al. (2018) with the exception of differences in nominotypical subspecies (discussed below) (Fig. 2). Similarly, difference in the alignments to Kadariya et al. (2018) resulted in different clusters in both networks. The sequence belonging to FM177759 was recovered in the smaller network between the Baluchistan and Nepalese black bears (Fig. 2), but its origin is unknown. Its high divergence to any other haplotype suggests that this could belong to a different subspecies. Two sequences (EF587265, KT964290) have ambiguous placement in the networks. The former is labeled as U. t. thibetanus in GenBank with unknown geographical origin, and the latter is a sequence from China, which grouped with Far East populations in Russia and Korea. The MCC tree and BI tree Fig. 2 Median-joining network of 54 Ursus thibetanus mtDNA control region sequences from mainland Asia and Taiwan with known locality of origin. Two U. t. japonicus haplotypes are included to show the position of this highly divergent clade relative to the other subspecies. The circles indicate individual haplotypes, and their sizes are proportional to the frequency of the haplotypes. Shared haplotypes (+) are as follows: cgrb2842 plus haplotypes from Russia (cgrb2056, cgrb2058-2060, cgrb838-842), Korea (cgrb54, crgb1728, cgrb2034-2036, cgrb2038-2041, Songwon43, Jangkang21, Duk-sung16, FJ895267, EF667005), and Russia/Korea (cgrb2055, cgrb672, cgrb843, cgrb2053); cgrb2037 plus a haplotype from Korea (Songwon9); NEP-A1 plus other haplotypes from Nepal (MG066704, MH281753, NEP-A2). Coding follows GenBank vouchers, when not available the GenBank accession number is used instead (Detailed information on haplotypes is given in Supplementary material, Table S1) were highly congruent in topology with high node posterior support probability (Pp) overall (Fig. 3). The Bayesian phylogenetic tree (Fig. 3) also suggested the existence of five clades/clusters, as had seen in the median-joining network analyses. The Baluchistan and Nepalese black bears were monophyletic (MCC, BI Pp = 1, ML = 74%) and are sister clade to all other subspecies from mainland and Taiwan. The GARLI ML 50% bootstrap consensus tree remained more unresolved with lower support overall.

Discussion
In this study, we sequenced a 338-bp fragment of the mitochondrial control region (CR) to assess the genetic diversity in Baluchistan black bear populations and a 680-bp to evaluate the phylogenetic relationships and divergence time between this population and other Asian populations. We admit that using a longer and more informative mtDNA sequences including full mitochondrial genome sequence can provide more robust estimates (Davison et al. 2011), but the dates proposed appear reasonable based on multiple sources of evidence (Waits et al. 1998).
Phylogenetic and median-joining network analyses of all U. thibetanus ssp. throughout mainland Asia (and Taiwan) recovered a highly diverged Baluchistan + Nepalese (Himalayan) Asiatic black bear population clade (Figs. 2, 3). Our findings suggest common Iranian and Nepalese ancestry at circa 120,000 years ago (Ky) (Fig. 3). These estimates interestingly coincide with major environmental changes in Eurasia following changes in climate from an interglacial period (characterized by the expansion of forests and animals, also known as Eemian or Riss-Würm interglacial), to a glaciation period (Hofreiterl and Stewart 2009). This interglacial period (130-115 Ky, MIS 5-2) was a period of warmer conditions that allowed forests to expand and replace tundra at higher latitudes, and opened new habitats to population expansions (Helmens 2013). Following the last glacial period, ice-cover throughout the Himalaya likely acted as a geographic barrier (Lei et al. 2014) that subdivided ancestral populations and led to genetic divergence. Consequently, this likely resulted in a constricted distribution of the Asiatic black bears into refugial areas, contributing to isolation of western populations that lived at the periphery of the subspecies range. Nevertheless, the dating range of the time since most recent common ancestor (TMRCA) of both populations falls within glacial and interglacial conditions (95% HPD: 234-33 Ky) and, therefore, accurate climatic conditions of such expansion (colder or warmer periods) cannot be postulated with confidence. Interestingly, fossil remains of the species (U. t. mediterraneus and U. t. permjak) have been found in the Caucasus and Europe dating to the middle Pleistocene (Baryshnikov 2010). These fossil records suggest that there have been several colonization events, most likely at interglacial periods towards the West.
The genetic diversity in Baluchistan black bears contrasts that of brown bears in Iran, which has been shown to be highly variable in terms of maternally inherited mtDNA (Ashrafzadeh et al. 2016). Unlike brown bears, the sampling recovered the lowest mitochondrial control region diversities  Table S1. Photo credit: F. Heydari reported in all the Asiatic black bears. Low mitochondrial genetic diversity in Baluchistan black bears may not be surprising as similar low-nucleotide diversity of mtDNA of the source population in the Himalayas has been previously reported (Kadariya et al. 2018), where only 3 haplotypes were recovered from 60 individuals. This can be partly explained by limited female-mediated gene flow as a result of female philopatry (female bears disperse and establish territories near their mothers) (Waits et al. 1998;Ohnishi and Osawa 2014) and/or by the small number of founders, as observed in the same species (Ishibashi and Saitoh 2004;Ohnishi et al. 2007Ohnishi et al. , 2009) and other bears (e.g., Zachos et al. 2008;Pérez et al. 2009;Swenson et al. 2011; but see Murphy et al. 2015Murphy et al. , 2018Murphy et al. , 2019. However, we admit that maternally inherited mtDNA markers may not be (and often are not) appropriate for making inferences about levels of genetic diversity/exchange at the population-level, because the species (as bears do in general) is characterized by male-biased dispersal (Waits et al. 1998). This is especially true when considering the fact that the dispersal of female bears to areas of marginal quality as those occupied by Baluchistan black bear habitat is very unlikely. An alternative explanation could be the decline of the population. Anthropogenic activities are known to have reduced Baluchistan black bear distributions and demographics (Yusefi 2013). Increasing sample sizes from less-sampled areas might challenge the current data for those areas, with the finding of additional haplotypes. Nevertheless, given the extremely small population of Baluchistan black bears in Iran (< 200 individuals), it likely represents most of the population genetic diversity. Furthermore, this sample size is also related to rareness of the species and sampling restrictions in remote and somehow unstable regions. It would be necessary to study the genetic diversity using other type of markers such as microsatellites, SNPs, etc., as such markers would allow us to better understand the genetic status of this endangered population.

Subspecies classification and taxonomy implications
Monophyly of Baluchistan + Nepalese populations questioning the validity of the U. t. gedrosianus subspecies described by Blanford (1877), calls for nominal taxonomic re-assessment and suggests that U. t. gedrosianus and the subspecies from central Nepal should be considered a single subspecies. A recent study (Kadirya et al. 2018) assigned the population from central Nepal to U. t. laniger (Pocock 1932), despite the fact that the correct name for Nepal populations should be U. t. thibetanus (Cuvier 1823) as it outlined by Pocock (1932). According to Pocock (1932), the range of the typical U. t. thibetanus includes Nepal (Terai: lowland belt that lies south of the outer foothills of the Himalayas) and extends to Myanmar and possibly southern China. Here, he also described U. t. laniger to the type locality in Kashmir, where Pakistan, India, and China meet and not Nepal or Bhutan as considered by Lan et al. (2017) and Kadirya et al. (2018). Indeed, the Asiatic black bear (U. thibetanus) was first found in the early nineteenth century in the mountains of Nepal (Cuvier 1823). However, Cuvier failed to assign a type locality for that taxon. Later, Pocock in 1932 allocated the typical thibetanus to Sylhet in Assam (northeastern India), because Cuvier gave the name thibetanus to a bear from that district. From this, we recognize thibetanus as the valid name for the Nepal population but not laniger. In the light of our results (i.e. monophyly of Baluchistan + Nepalese populations), and depending on further studies using autosomal nuclear markers (microsatellites and SNPs), and morphology, these two black bear populations might be considered as synonymous under the prior U. t. thibetanus trinomial (with gedrosianus just as junior synonym) as prescribed by the rules of the International Code of Zoological Nomenclature (ICZN), (ICZN 1999). Last, it is unclear as to which subspecies the population from Pakistan (Lan et al. 2017) belongs to; however, the preliminary analysis of the short control region (180 bp) fraction was identical to the Iranian haplotype. Most importantly, our study supports new evidence for only five subspecies (thibetanus, formosanus, japonicus, mupinensis, and ussuricus) rather than seven, as suggested by previous works (e.g. IUCN; Wozencraft 2005).

Conservation implications
Despite once considered locally extinct and rediscovered in the 1970s, the Baluchistan black bear has received little attention from the Iranian Department of Environment (DoE), a legal entity responsible for issuing the national directives for environment and wildlife protection. Limited surveys have been carried out for the species and most of its range remains unprotected. Conservation for this endangered species is rarely enforced, even though this species is currently classified as "protected" by national Iranian law (Yusefi et al. 2010;Yusefi 2013). Worryingly, over the past decades the distribution of this bear has been reduced, becoming locally extinct in some areas (e.g., Khabr National Park and Birk Protected Area) (Ziaie 2008). Furthermore, because of the scarcity of suitable habitat, distribution patchiness, and geographically separated subpopulations, it is very likely that these areas hold only a few individuals. Such fragile subpopulations are at risk due to both genetic impoverishment and demographic threats. Immediate conservation actions such as reducing human-bear conflict, designation of new protected areas and protection of corridor habitats may help these small subpopulations survive.
Even with limited sampling, our study sheds light on the phylogenetic affinity of Baluchistan black bears, and (to some extent) on the taxonomy of Asiatic black bears below the species level, as well as providing a genetic benchmark for future molecular analyses of one of the least-known populations of the black bears in Asia. Nevertheless, we do recognize that our results should be considered with caution due to the limitations of our study.