Genetic diversity within and among gelada (Theropithecus gelada) populations based on mitochondrial DNA analysis

Abstract

We studied genetic characteristics within and among gelada (Theropithecus gelada) populations inhabiting the southern and northern plateaus of Ethiopia. Twenty-one mtDNA haplotypes were identified. Geladas on the southern plateau were genetically separated from those on the northern plateau, with a large differentiation as indicated by F_st values of 0.665–0.917. The difference between the subspecies (T. g. gelada and T. g. obscurus) on the northern plateau highlighted a substantial genetic variation. Divergence times were estimated as ~250000 years between subspecies on the northern plateau and ~400000 years between those on the northern and southern plateaus. The genetic differentiation between the geographically distant Simien Mountain and Arsi groups (640 km) was ~2/3 of that between the geographically closer Debra Libanos and Arsi groups (250 km). The difference between subspecies within the northern plateau was similar to that between Papio hamadryas hamadryas and P. h. anubis. The difference between geladas in the north and south was similar to inter-subspecies differences in other mammals, and thus suggested that the Arsi geladas belong to subspecies T. g. arsi (tentative name). Considering the present distribution, geladas appear to have a complicated history of speciation. However, further analyses based on genetics, morphology and ecology are required to confirm these findings.

Introduction

The gelada (Theropithecus gelada) is a primate species indigenous to Ethiopia that was first discovered on the northern Ethiopian plateau and was described by Rüppell (1835). A subspecies was later documented by Heuglin (1863). Geladas from the northern (north of Lake Tana and west of the Takazze river) and southern regions (south of Lake Tana and east of the Takazze river) are believed to be T. gelada gelada and T. gelada obscurus, respectively; however, a clear distinction between them has not yet been determined (Yalden et al., 1977, 1996; Gippoliti, 2010). Several studies have detailed the natural history of this species (Jolly, 1972; Delson, 1993; Jablonski, 1993). Recently, Frost et al. (2014) found fossils of early Theropithecus from in the Afar area, dating to 3.6–3.8 million years ago (mya). Many studies regarding the social and ecological aspects of geladas have been conducted since the late 1970s (Dunbar and Dunbar, 1975; Kawai, 1979; Bergman and Beehner, 2008; Bergman et al., 2009) mainly in the Simien Mountains on the northern plateau. Furthermore, Snyder-Mackler et al. (2014) studied in detail the individual relationship among gelada populations with stratified social structure via a genetic approach. Near Debre Libanos, the southernmost region of the northern plateau, Shotake and Nozawa (1984) conducted a genetic study of variation in blood proteins as markers. On the basis of low genetic variations in four neighboring groups near Debra Libanos, they demonstrated that gelada species had been subjected to a bottleneck effect during speciation. Matsubayashi et al. (1998) measured blood erythropoietin levels and red blood cell counts in geladas from the Simien Mountains and reported that they were 1.5- to 2-fold higher than those of geladas from the plains. They concluded that these alterations were caused by physiological adaptations to the high altitude rather than by genetic events.

Mori and Belay (1990) discovered and described a new group of geladas in the Arsi region, which is bordered by the uppermost stream of the Wabe-Shebelle river on the southern plateau. However, Groves (1972) had previously assumed that some geladas were present in this area and named this virtual subspecies senex (Yalden et al., 1977; Mittermeier et al., 2013). Gippoliti (2010) pointed out that the subspecies name ‘senex’ is nonsense because it is already a junior synonym of gelada. Iwamoto et al. (1996) and Mori et al. (1997, 1999, 2003) later conducted detailed studies on the social and ecological behaviors of this group and documented differences in social structures, behaviors, and fur coloration, as compared with geladas on the northern plateau. Specifically, Arsi gelada single-male units and bands are characterized by much smaller numbers than those of northern plateau geladas, and their behaviors are characterized by infanticide, tree climbing, and mobbing behavior in response to predators, which are not observed on the northern plateau. Furthermore, their fur color is lighter. Belay and Shotake (1998) compared genetic variation in blood proteins of these geladas with that of those from the northern plateau and reported large genetic differentiation, estimated to correspond to about 0.4 million years of divergence time by Nei’s genetic distance (Nei, 1972) between the groups, and suggested the third subspecies. Belay and Mori (2006) suggested, on the basis of a preliminary mitochondrial DNA (mtDNA) analysis, that the Arsi gelada is a distinct subspecies. Gippoliti (2010) described ambiguities in the current classifications and distributions of geladas and emphasized that detailed surveys and studies are required for conservation of geladas.

Based on the results of these previous studies, in the present study, we analyzed mtDNA to identify genetic differences between the northern and southern groups of geladas and to explain the regional subspeciation within the northern plateau. We attempted to establish the phylogeographical features of this species in view of the subspeciation by associating the genetic characteristics with estimated times of divergence.

Materials and Methods

Blood sample collection and mtDNA analysis

Blood was collected from a total of 113 individuals from the following seven populations from 1978 to 1997: Simien Mountain (NSi), Magdera (NMa), Goshmeda (NGo), Debra Libanos (NDb), Shinkrto (NSh), Fiche (NFi), and Arsi (SAr), where the abbreviations N and S in parentheses indicate the northern and southern plateaus (Table 1, Figure 1). The guidelines for field research of the Kyoto University Primate Research Institute and the rules and regulations of the Ethiopian government at that time were followed during blood collection: we trapped animals in nets and anesthetized them with ketamine chloride doses of 1 ml/5 kg. From each animal, 1–5 ml of blood was withdrawn from a vein in the arm, and the animal’s sex and age were determined. Animals were released after sedation had completely worn off. For mtDNA analysis, the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) was used to extract genomic DNA from frozen blood samples collected by the author (T.S.) during 1978–1997 for blood protein analysis. The 390 bases (including one deletion) corresponding to positions 15462–15851 of D-loop region of Papio hamadryas (GenBank accession number Y18001) were sequenced. In T. gelada, this region corresponds to 389 bases at positions 15467–15855 (GenBank accession number FJ785426). The reaction mixture for polymerase chain reaction (PCR) contained genomic DNA, 10 pmol of each primer, 2 μl of 2.5 mM dNTPs, 1 U of Ex Taq polymerase, and 2.5 μl of 10× Ex Taq buffer (TaKaRaBio, Inc., Tokyo, Japan). The thermocycling conditions were as follows: initial denaturation at 94°C for 2 min, followed by 35 cycles at 94°C for 30 s, annealing at 50°C for 35 s, and extension at 72°C for 2 min, and a final extension at 72°C for 4 min. Two primers, outF (forward) 5′-CTGGCGTTCTAACTTAAACT-3′ and outR (reverse) 5′-GTAGTATTACCCGAGCGG-3′, were used to obtain a PCR fragment spanning the D-loop of the mtDNA genome (Hapke et al., 2001). PCR products were purified using the Gene Clean Kit (Funakoshi Co., Ltd., Tokyo, Japan), and the size and approximate yield were determined by standard agarose gel electrophoresis. The purified PCR products were cycle-sequenced using the BigDye Terminator v. 3.1 Cycle Sequencing Kit and the ABI PRISM^® 377 DNA Sequencer (Applied Biosystems, Waltham, MA, USA).

Table 1 Sampling populations and sites of geladas in this study

Population with abbreviation	Nos. of animals sampled	Approximate number of animals in original population	Position	Sampling year
Simien mountain (Nsi)	21	412	13.1655N, 38.0726E	1995
Near Magdera (NMa)	10	210	11.1121N, 39.1412E	1997
Goshmeda (NGo)	9	413	9.4526N, 39.4320E	1993
Debra Libanos (NDb)	14	168	9.4301N, 38.5006E	1978
Shinkrt (NSh)	10	204	9.4406N, 38.4851E	1978
Fiche (NFi)	10	201	9.4934N, 38.4341E	1979
Arsi (SAr)	39	169	7.5033N, 39.5130E	1996

Figure 1

Populations and sites sampled in this study with reference to maps (scale 1/50000 or 1/500000) published by the Ethiopian mapping agency in 1985 and Google Earth.

Data analyses

All D-loop sequences of geladas from the seven populations examined in this study were deposited in the DNA Data Bank of Japan, European Molecular Biology Laboratory database, and GenBank database under accession numbers LC018113–LC018133. The sequences were multiply aligned with ClustalW (Larkin et al., 2007). Nucleotide diversity (Nei, 1987) and F_st (differentiation index between populations) were calculated using the Arlequin version 3.0 (Excoffier et al., 2005). A maximum-parsimony haplotype network was generated with Network 4.6.1.2 Phylogenetic Network Software based on a median-joining network (Bandelt et al., 1999). A phylogenetic tree was constructed using Kimura’s two-parameter distance (Kimura, 1980) with MEGA 5.1 software (Tamura et al., 2011) by the neighbor-joining (NJ) method (Saitou and Nei, 1987). To test past population expansion, Fs statistics (Fu, 1997), Tajima’s D (Tajima, 1989), and the raggedness index (R) (Harpending, 1994) were calculated using Arlequin, as well as by estimation of demographic parameters for the distribution of pair-wise differences (mismatch distribution) (Rogers and Harpending, 1992; Schneider and Excoffier, 1999).

Estimation of divergence time

Divergence times of the mtDNA haplotype groups were estimated using the BEAST 2.1.3 software package for Bayesian evolutionary analysis (Drummond and Bouckaert, 2015). As calibrations we used two divergence points in the analysis. One is the fossil data between the divergence of the genera Theropitecus and Papio (assumed to be ~5.00 ± 1.00 mya; Frost et al., 2014) and the other is blood protein data between the divergence of geladas from the northern and southern plateaus (~0.40 ± 0.05 mya; Belay and Shotake, 1998). In the time estimation, MCMC and burn-in conditions were set to 10000000 runs and 1%, respectively. Finally, we constructed the consensus tree and DensiTree (Bouckaert, 2010), depicting mtDNA evolution among the study populations.

Results

mtDNA analysis

Distribution of haplotypes

Table 2 shows the nucleotide sequences of 21 haplotypes that were polymorphic at 75 sites. H1–H6 and H7–H19 were observed in the Simien Mountain (NSi) and Blue Nile basin system populations (NMa, NGo, NFi, NSh, and NDb), respectively, whereas H20–H21 were observed in the population from the southern plateau of Arsi (SAr) across the Rift Valley.

Table 2

Observed 21 haplotypes, the numbers of each haplotype and populations and their 389 nucleotide aliment positions according to the D-loop of Theropithecus gelada (FJ785426; 15467~15855). 75 polymorphic sites exist among 04 (15470) to 360 (15826).

The nucleotide sequences of the haplotypes in Table 2 show that each haplotype was unique to one of the three regions. Transversions were observed at positions 73 (A-C) and 305 (C-T and A). Figure 2 depicts the network constructed on the basis of this table. The result showed a clear division into three segments. Individuals indicated by ♂ (H5, 6, 8, 14, 17, 18, 19 and 20) are young or adult males. H5, H17 and H18 are, in particular, greatly differentiated from the mother groups. These results suggest that males migrate between populations of geladas, as occurs in other papionins. The composition list of each population is shown in the Appendix.

Figure 2

Minimum spanning haplotype network constructed using Network 4.6.1.0 based on the median-joining algorithm. The symbol ♂ indicates adult or young males. The sizes of circles are proportional to the numbers of haplotypes.

The NDb and NSh populations shared three common haplotypes. Despite the close physical proximity (10 km) between the NFi and NSh populations, there were no haplotypes in common, possibly a consequence of the small number of samples collected from these groups. The NMa population, which is located in the northernmost region among the populations inhabiting the Blue Nile basin system, is greatly differentiated from the other populations. The NGo population also has two unique haplotypes that are slightly diverged from the other populations. Although 39 geladas in the SAr population were analyzed, only two haplotypes with one substitution were found. The SAr population on the southern plateau was more closely related to the NSi population in the Simien Mountains, which is located in the northernmost regions, by 30 substitutions (including one transversion) than the NMa, NGo, NFi, NSh, and NDb populations of the Blue Nile basin region (Figure 2). However, H21 of the SAr population had 44 nucleotide substitutions (including two transversions) relative to H11 in the Blue Nile basin (Table 2). Similar trends were observed in several parameters to be described later. Figure 3 shows a schematic representation of the NJ tree of haplotypes, with P. hamadryas (GenBank accession number Y18001) represented as an outgroup. This figure also shows a clear division of the haplotypes into three regions.

Figure 3

NJ phylogenetic tree of 21 gelada haplotypes with an out-group sequence of Papio hamadryas (GenBank accession number Y18001) constructed using MEGA software. Number of nucleotide substitutions per site as indicated by the scale are the Kimura’s two-parameter distances. The number of each branch indicates the bootstrap value (>50) after 1000 replications.

Statistical analyses of haplotype variations

Next, several genetic parameters were calculated from Table 2 and are shown in Table 3. Except for the NDb and NSh populations, each group had unique haplotypes (Uh). The NDb and NSh populations have both unique and common types. There was a high degree of haplotype and nucleotide diversity in the NSi, NDb, and NSh populations. The NDb and NSh populations are neighboring groups presumed to have diverged from a single population (Shotake and Nozawa, 1984). Table 4 summarizes the pairwise F_st (differentiation index) between populations. The highest F_st values (0.665–0.917) were observed for the SAr population, which was located across the Rift Valley compared with those of the northern plateau populations, suggesting long-term isolation. In Table 4, the smallest F_st value of the SAr population was 0.665 for the NSi population in comparison with those of the other six populations. Arsi geladas were more genetically similar to Simien than Blue Nile geladas. Moreover, no genetic differences between the NDb and NSh populations were observed.

Table 3 Summary statistics of population diversity indices in seven populations of gelada.

Population	n	P	H	Uh	Hd ± SD	π ± SD
NSi	21	28	6	6	0.767 ± 0.069	0.0188 ± 0.0102
NMa	10	11	2	2	0.200 ± 0.154	0.0057 ± 0.0039
NGo	9	3	2	2	0.500 ± 0.130	0.0039 ± 0.0009
NSh	10	18	5	2	0.755 ± 0.130	0.0211 ± 0.0038
NDb	14	20	4	1	0.747 ± 0.066	0.0221 ± 0.0023
NFi	10	2	3	3	0.600 ± 0.131	0.0017 ± 0.0016
SAr	39	1	2	2	0.051 ± 0.048	0.0001 ± 0.0001
All populations	113	75	21	18	0.858 ± 0.026	0.0691 ± 0.0338

n, number of samples; P, number of polymorphic sites; H, number of haplotypes; Uh, unique haplotype; Hd, haplotype diversity; π, nucleotide diversity.

Table 4 Genetic differentiation (F_st) indicated among seven populations of gelada.

Population	NSi	NMa	NGo	NSh	NDb	NFi	SAr
NSi	—
NMa	0.456	—
NGo	0.339	0.656	—
NSh	0.235	0.522	0.368	—
NDb	0.239	0.497	0.362	0.000*	—
NFi	0.301	0.600	0.448	0.322	0.321	—
SAr	0.665	0.917	0.852	0.729	0.770	0.816	—

All comparison values were significant at P < 0.0001, except ^* P = 0.621.

We conducted molecular variance analysis of the gelada populations (Table 5). Based on these data analyses, the entire gelada population was divided into three major regions: the Simien Mountain, Blue Nile basin, and Arsi regions. Accordingly, further analysis was performed and showed that most (77.1%) of the variance was attributable to the diversity among these three regions, in strong agreement with the results of the previous analysis. Some variation (13.98%) was found within the Blue Nile basin and Simien Mountain regions.

Table 5 Analysis of molecular variance (AMOVA) of seven populations of gelada classified into three regions defined by geographical locality

Source of variation	df	Sum of squares	Variance component	Percent of variation
Among regions	2	1196.65	14.99	77.10
Among population within regions	4	121.45	2.72	13.08
Within populations	106	183.99	1.78	8.92

Note: Region 1, NSi; Region 2, NMa, NGo, NFi, NSh and NDb; Region 3, SAr.

Population expansion

Differences in phylogeographic distribution of mtDNA haplotypes were considerable within and among the three regional gelada populations. Tests for population expansion with Fu’s measure and Tajima’s D of selective neutrality were not significant for any regional population. In mismatch distribution and raggedness analyses, Simien Mountain and Blue Nile basin populations did not show significant values, whereas the value for the Arsi population seemed to be significant (Table 6). However, the Arsi population showed only one substitution and we judged this differentiation to be negligible. Overall, we could not from these estimates identify population expansion in any regional population.

Table 6 Summary statistics for the population expansion test

Populations	Simien	Blue Nile	Arsi
Sample size	21	53	39
No. of substitution sites	28	31	1
Mismatch distribution parameter
τ	10.313	15.996	3.000
θ0	0.000	0.000	0.053
θ1	10.596	17.276	0.053
Raggedness index (R)	0.171	0.034	0.808
P (Sim.R ≥ Obs.R)	0.000	0.090	0.860
Tajima’s D	−0.237	1.034	−1.126
P (D sim < D obs)	0.457	0.885	0.143
Fu’s test (Fs)	4.703	3.166	1.429
P (sim_Fs < obs_Fs)	0.965	0.861	0.042

Note: τ, the time after expansion in mutational units; θ₀, the mutation parameter before expansion; θ₁, the mutation parameter after expansion. For details of parameters, see Rogers and Harpending (1992) and Schneider and Excoffier (1999).

Divergence time

The divergence time among three regions of geladas (Simien Mountain, Blue Nile basin system, and Arsi) was estimated and is shown in Figure 4. The divergence time of the most recent common ancestor (MRCA) of geladas was estimated as 0.4057 ± 0.0934 mya (node 2), and the divergence time of the MRCA of the northern plateau geladas was estimated as 0.2474 ± 0.1327 mya (node 3). Haplotype groups among the Blue Nile basin were divided into two groups ~0.098 mya (node 5), and this estimate can be attributed mainly to the NMa population, which is isolated by a geographical barrier. Figure 5 shows some connections between the Simien Mountain and Arsi populations described above.

Figure 4

Consensus tree of the haplotypes of gelada populations with Papio hamadryas as the outer group. Gray bars indicate 95% confidence intervals of divergence time.

Figure 5

DesiTree (all + root canal) showing the posterior probability of 10000 trees illustrated by Bayesian divergence time analysis using Beast 2. Trees are depicted in green and the root canal is shown in blue.

Discussion

Relationship between haplotype and geographical distribution

In Table 3, the high haplotype diversity in the NSi population is likely due to the large size of this population (Beehner et al., 2007; Snyder-Mackler et al., 2014). Snyder-Mackler et al. (2014) researched in detail the kin-relationships within gelada populations based on mtDNA sequenced haplotypes as markers around the Sankaber area in the Simien Mountains. They sequenced 409 bases of D-loop which corresponded to nucleotide positions 15438–15846 of the gelada mitochondrial genome (GenBank accession number FJ785426) and recognized 15 haplotypes, of which six haplotypes were configured from six animals. Our sequence corresponded to nucleotide positions 15467–15855 of the mitochondrial genome and overlapped 380 bases with those sequences. All polymorphic sites of their haplotypes are distributed within the overlapping sequence zone (380 bases). Our haplotypes observed in the NSi population except H2 were included in their haplotypes. Thus, the genetic difference of populations between the Gich plateau, our study site, and around the Sankaber area, which is 10–25 km along the cliff, is small as compared with the relation of our NDb and NSh populations. We can estimate that a large breeding population of gelada spread over the Simien Mountains, where there are many connecting cliffs with gelada habitats.

As suggested by the data in Table 3, the high diversities observed in the NDb and NSh populations are also reasonable if we consider that the populations are closely connected with a large reproductive group of several thousand geladas that seems to inhabit both sides of the basin of the Jemma and Wonchit rivers. Kifle et al. (2013) counted 1465 geladas in only one location in the Wonchit valley (10.15N, 39.17E) shown as W. valley in Figure 1. Also, in 1991 we observed approximately 150 geladas in the Mendida valley (9.39N, 39.23E) shown as M. Despite the close proximity between the NFi and NSh groups, they had no haplotypes in common, possibly owing to the small number of samples collected. Shotake and Nozawa (1984) reported that the NFi and NSh populations share a few variants of blood protein genes. This suggests that there are some immigrants between both populations. The presence of haplotypes next to each other (H7 and H12) with only two substitutions could represent a recent isolation event.

The NMa population has greatly diverged from the other populations. The presence of a male with 11 substitutions (H18) suggests that geladas around this region were distributed over a large area, as was the case for the NSi population. Yalden et al. (1977) described that there were many gelada habitats along the Bishlo river, including where the NMa population is currently located. In the southwest area of Dese is a chain of mountains with four aligned peaks that exceed 3800 m, stretching from Mount Amba Farit (4247 m) toward the east and ending at the cliff of the Rift Valley. This chain is a watershed of the Bishlo, Wonchit, and Weiqua rivers but has no cliffs. This region is probably the blank area of the distribution that was indicated by Yalden and Largen (1992) and Gippoliti (2010). This geological condition may be an isolation factor that contributed to the subdivision of the populations of the Jemma or Wonchit rivers located to the south.

The NGo population also shows two unique haplotypes, indicating little divergence from the other populations that may be attributable to a geographical factor in the area that they inhabit, where the branch of the Jemma river originates. This region, which is located downstream of the river, has no cliffs and is surrounded by cropland, so that recent migration is impossible. Considering that this group shows little divergence from downstream groups (H10 and H11 to H7, H12, and H13 in Figure 2) and geographical factors, the differentiation event is believed to have occurred recently.

The cliffs of the Rift Valley begin 50 km northeast of Addis Ababa and continue to the north for 600 km. Several populations inhabit this area and form a large reproductive group. Thirty kilometers north of Guassa (10.18N, 39.48E in Figure 1) occurs the origin of a branch of the Wonchit river, which forms cliffs that are inhabited by many geladas (Fashing et al., 2010; Ashenafi et al., 2012). This area is adjacent to the Rift Valley, which is accessible from downstream. This geographical factor is believed to have contributed to the diversity of the populations around the Blue Nile basin. The SAr population showed extremely low genetic diversity, as mentioned above. Similarly to the findings of Belay and Shotake (1998) about genetic variation in blood proteins, our results support the assumption that this group lost variation in a bottleneck, shrank, and then developed into the present group after migration from the northern plateau.

Recognition of subspecies

The classification of T. gelada gelada and T. g. obscurus is under debate (Yalden and Largen, 1992; Yalden et al., 1996; Gippoliti, 2010). According to our observations, the boundary between these subspecies is the central tableland, crossed by the road from Wereta to Weldya, on the northern plateau (a watershed of the Takazze and Blue Nile rivers, with Mount Guna, 4275 m high, on the west). The Takazze and branches of the Blue Nile intertwine at the center of this plateau. Although parts of these branches shaped the mountain areas without cliffs, most of the regions are plateaus with elevations of 2500–3000 m. Two-thirds of the area, west of the watershed, was presumed to be the border between the two subspecies. The east end of this plateau continues to the cliffs in the Rift Valley, which stretches into the north of the plateau. In this region, a southeast branch of the Takazze river originates, and the geladas inhabit the cliffs of this region. Mount Abune Yosef (4190 m) is located in the northeast, and the geladas inhabit the northwest cliffs of this mountain massif. According to our observations, the gelada from this region have dark fur and a larger triangular pink chest patch surrounded by clear white fur, resembling the appearance of those of the Blue Nile basin and matching those from the obscurus habitat (east of the Takazze river). However, because of the large size of the area where the samples were collected, we cannot simply conclude that this water basin system has a greater diversity than that of other areas. Further analyses of the other groups should be performed in future studies.

The results of this study showed significant differentiation between the gelada and obscurus subspecies. Figure 2 and Table 2 show that a minimum of 27 substitutions including one transversion were detected between the haplotypes of the NSi population, which was considered to represent T. g. geladas, and those of the NFi or NGo populations, which were considered to represent T. g. obscurus. Furthermore, a minimum of 30 substitutions including one transversion were detected between H2 from the NSi population and H21 from the SAr population on the southern plateau. These populations were the furthest apart geographically and were separated by the Rift Valley. The differentiation index F_st was 0.235 for the NSi and NSh populations, which were closest to each other, and 0.665 for the NSi and SAr populations. Northern plateau populations and the SAr population showed extremely high F_st values (0.665–0.917). If the gelada and obscurus are subspecies (Gippoliti, 2010; Mittermeier et al., 2013), then the highly differentiated populations on the southern and northern plateaus should also be treated as subspecies T. g. arsi, which we named tentatively. Hapke et al. (2001) estimated an F_st of 0.599 between P. hamadryas hamadryas and P. h. anubis in the same region (a 342 base sequence), in which we examined 389 bases of the D-loop of mtDNA. Thus, the genetic difference (F_st = 0.665) between the southern and northern plateau populations was at the same level as the differentiation between P. h. hamadryas and P. h. anubis. In addition, Shotake (1981) estimated a divergence time of approximately 0.35 mya between P. hamadryas and P. anubis by blood protein markers. This value is comparable with that between Arsi and northern plateau geladas (Belay and Shotake, 1998).

Moreover, with respect to the differentiation of the southern and northern plateaus populations, the SAr group appeared to be more closely related to the NSi group than to populations inhabiting the Blue Nile basin, which is closer in geographical distance. Furthermore, the DensiTree in Figure 5 suggests that there may have been a connection between the Arsi and Simien Mountain regions.

Mori and Belay (1990) noted that the external features of Arsi geladas were closer to the description by Yalden et al. (1977) of Simien gelada (gelada) rather than that of Debra Sina gelada (obscurus), which is a neighboring troop of our NGo population. We accordingly reviewed the photographs from past expeditions. We propose that the representative fur colors from each region, which suggested that the fur coloration of the SAr group was more similar to that of the NSi group than to those of the NDb and NFi groups. However, it will be necessary to quantify the color. Arsi geladas have a lighter overall color than the NSi group. For this reason, the name obscurus was assigned to this subspecies rather than gelada, which likely means ‘blackish.’ Gippoliti (2010) also mentioned that ‘Southern gelada’ (T. g. obscurus) were similar to ‘Northern gelada’ (T. g. gelada), but members of this group have a darker mane and extensive, pure white fur surrounding the pink chest patch.

Divergence time

Given that there may be some problems associated with estimating the differentiation time and environmental factors between populations using only mtDNA markers (Zachos et al., 2013), we also assessed nuclear DNA (blood proteins) to calibrate time. We presented one tentative example of the relationship between the differentiation of gelada populations and environmental factors.

The Great Rift Valley was formed more than 2.5–3.0 mya (Baker et al., 1972). If this valley is a factor in the isolation between the northern and southern populations of gelada, our estimate of 0.4057 mya appears too small. However, the Rift Valley did not suddenly appear, but rather was gradually formed by crustal movement and volcanic activity. The geladas may have migrated from the northern plateau to the southern plateau after formation of the Rift Valley 0.4 mya (Iwamoto et al., 1996). After these events, for example, increasing volcanic activity may have occurred, as evidenced by the many small fresh craters between Nazaret and Lake Bassaka in Figure 1. This area is the narrowest part of the Rift Valley bottom, which appears to be a potential migration route.

The divergence time between T. g. gelada and T. g. obscurus was estimated as 0.247 mya (Figure 4). We did not find evidence of a drastic crustal movement around this boundary. However, several glacial periods in the area expanded the upper streams of the Takazze and Blue Nile rivers, and changes in the flora and fauna indicate these transitions (Hamilton, 1982). This area receives considerable rainfall (~1500 mm: Wolde, 1970; Gamachu, 1977) during the 3 month rainy season. We propose that these heavy rains have continued to carve the earth and formed the large branches of both rivers, leading to the loss of all cliffs. Thus, these events may be the factors responsible for the isolation of the two subspecies. Haplotypes along the Blue Nile basin may have been divided into two groups ~0.098 mya and this estimate can be explained by the geographical barrier without cliffs between the NMa group and other populations mentioned earlier. Finally, based on the above genetic results and those of our predecessors (Mori and Belay, 1990; Iwamoto et al., 1996; Mori et al., 1997, 1999, 2003; Belay and Shotake, 1998; Belay and Mori, 2006), we propose that Theropithecus gelada should be classified into three subspecies: gelada, obscurus, and arsi, the last of which is named tentatively. However, considering that species ranks are assigned to P. hamadryas and P. anubis in the latest comprehensive review of primate taxonomy (Mittermeier et al., 2013), a species rank might be more appropriate for arsi instead of the subspecies rank as suggested by Gippoliti (2010).

Acknowledgments

We thank Profs. A. Mori, T. Iwamoto, M. Kawai, K. Nozawa, and K. Matsubayashi as well as Drs. G. Belay, M. Shimada O. Hishida, A. Yamane, and S. Yrga for their support of this long-term research project. We thank Dr. M. Huffman for revising English of some parts in the manuscript. We would also like to thank the staff members of Addis Ababa University, particularly Profs. E. Bekele, Y. Mekonen, A. Bekelle, and B. Petros, and Mr. L. Abune of the Ethiopian Wild Life Conservation Organization for their assistance in this joint study. We also thank Mrs. S. Kawamoto and R. Ebihara for technical assistance. This work was supported by Grants-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan (Nos. 02041104, 04041084, 06041065, 15570194) and 1978, 1983, 1999, and 2002 Grant of Overseas Special Research Program of Primate Research Institute, Kyoto University from the Ministry of Education, Culture, Sports, Science and Technology of Japan. Finally, we deeply thank three anonymous reviewers and the Editor-in-Chief for fruitful and constructive comments on this paper.

Ethical standard

This work is a continuous joint study with the Ethiopian Wildlife Conservation Organization and Addis Ababa University since 1973. Research methods were approved by the Guideline for Animal Care and Use Committee of Kyoto University and the procedures of blood collection were performed according to regulations of the Ethiopian government at that time.

Appendix

Appendix Investigated populations and animals’ ages, sexes, and haplotypes

Simien Mountain (NSi)
NSi 01	Am	haplo 6
NSi 02	Am	haplo 2
NSi 03	Am	haplo 4
NSi 04	Am	haplo 4
NSi 05	Am	haplo 1
NSi 06	Af	haplo 1
NSi 07	Af	haplo 2
NSi 08	Af	haplo 1
NSi 09	Af	haplo 1
NSi 10	Af	haplo 1
NSi 11	Am	haplo 5
NSi 12	Ym	haplo 1
NSi 13	Yf	haplo 3
NSi 14	Jm	haplo 1
NSi 15	Ym	haplo 2
NSi 16	Ym	haplo 3
NSi 17	Jm	haplo 4
NSi 18	Jm	haplo 1
NSi 19	Ym	haplo 1
NSi 20	Jm	haplo 3
NSi 21	Jf	haplo 4

Magdera (NMa)
NMa 01	Am	haplo15
NMa 02	Am	haplo18
NMa 03	Am	haplo15
NMa 04	Af	haplo15
NMa 05	Af	haplo15
NMa 06	Af	haplo15
NMa 07	Ym	haplo15
NMa 08	Ym	haplo15
NMa 09	Ym	haplo15
NMa 10	Af	haplo15

Goshmeda (NGo)
NGo 01	Ym	haplo10
NGo 02	Ym	haplo11
NGo 03	Yf	haplo10
NGo 04	Yf	haplo10
NGo 05	Yf	haplo11
NGo 06	Yf	haplo10
NGo 07	Inm	haplo10
NGo 08	Inf	haplo10
NGo 09	Bam	haplo11

Fiche (NFi)
NFi 01	Af	haplo 7
NFi 02	Jf	haplo 7
NFi 03	Jf	haplo 7
NFi 04	Jf	haplo 9
NFi 05	Jm	haplo 7
NFi 06	Ym	haplo 9
NFi 07	Ym	haplo 8
NFi 08	Ym	haplo 7
NFi 09	Af	haplo 7
NFi 10	Jm	haplo 9

Shinkrto (NSh)
NSh 01	Am	haplo16
NSh 02	Ym	haplo16
NSh 03	Jf	haplo16
NSh 04	Jf	haplo16
NSh 05	Ym	haplo19
NSh 06	Ym	haplo12
NSh 07	Jf	haplo13
NSh 08	Af	haplo13
NSh 09	Ym	haplo14
NSh 10	Yf	haplo16

Debra Libanos (NDb)
NDb 01	Af	haplo16
NDb 02	Am	haplo13
NDb 03	Am	haplo13
NDb 04	Ym	haplo12
NDb 05	Af	haplo13
NDb 06	Ym	haplo16
NDb 07	Ym	haplo17
NDb 08	Af	haplo13
NDb 09	Yf	haplo12
NDb 10	Yf	haplo16
NDb 11	Ym	haplo16
NDb 12	Jf	haplo12
NDb 13	Ym	haplo16
NDb 14	Jm	haplo13

Arsi (SAr)
SAr 01	Am	haplo 20
SAr 02	Am	haplo 21
SAr 03	Am	haplo 21
SAr 04	Am	haplo 21
SAr 05	Am	haplo 21
SAr 06	Jf	haplo 21
SAr 07	Ym	haplo 21
SAr 08	Ym	haplo 21
SAr 09	Ym	haplo 21
SAr 10	Am	haplo 21
SAr 11	Am	haplo 21
SAr 12	Ym	haplo 21
SAr 13	Am	haplo 21
SAr 14	Ym	haplo 21
SAr 15	Ym	haplo 21
SAr 16	Ym	haplo 21
SAr 17	Ym	haplo 21
SAr 18	Ym	haplo 21
SAr 19	Am	haplo 21
SAr 20	Ym	haplo 21
SAr 21	Jm	haplo 21
SAr 22	Jm	haplo 21
SAr 23	Yf	haplo 21
SAr 24	Am	haplo 21
SAr 25	Ym	haplo 21
SAr 26	Ym	haplo 21
SAr 27	Am	haplo 21
SAr 28	Am	haplo 21
SAr 29	Am	haplo 21
SAr 30	Ym	haplo 21
SAr 31	Jf	haplo 21
SAr 32	Jf	haplo 21
SAr 33	Af	haplo 21
SAr 34	Af	haplo 21
SAr 35	Af	haplo 21
SAr 36	Af	haplo 21
SAr 37	Af	haplo 21
SAr 38	Ym	haplo 21
SAr 39	Am	haplo 21

Note

Af, adult female

Am, adult male

Yf, young female

Ym, young male

Jf, juvenile female

Jm, juvenile male

Inf, infant female

Inm, infant male

Bam, baby male

References

•  Ashenafi  Z.T.,  Leader-Williams  N.L., and  Coulson  T. (2012) Consequences of human land use for an Afro-alpine ecological community in Ethiopia. Conservation and Society, 10: 209– 216.
•  Baker  B.H.,  Mohr  P.A., and  Williams  L.A.J. (1972) Geology of the Easter Rift System of Africa. Geological Society of America Special Papers, 136: 1–68.
•  Bandelt  H.J.,  Forster  P., and  Rohl  A. (1999) Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution, 16: 37–48.
•  Beehner  J.C.,  Berhanu  G.,  Bergman  T.J., and  Cann  M.C. (2007) Population estimate for geladas (Theropithecus gelada) living in and around the Simien Mountains National Park, Ethiopia. Ethiopian Journal of Science, 30: 149–154.
•  Belay  G. and  Mori  A. (2006) Intraspecific phylogeographic mitochondrial DNA (D-loop) variation of gelada baboon, Theropithecus gelada, in Ethiopia. Biochemical Systematics and Ecology, 34: 554–561.
•  Belay  G. and  Shotake  T. (1998) Blood protein variation of a new population of gelada baboon (Theropithecus gelada), in the southern Rift Valley, Arsi Region, Ethiopia. Primates, 39: 183–198.
•  Bergman  T.J. and  Beehner  J.C. (2008) A simple method for measuring color in wild animals: validation and use on chest patch color in geladas (Theropithecus gelada). Biological Journal of the Linnean Society, 94: 231–240.
•  Bergman  T.J.,  Ho  L., and  Beehner  J.C. (2009) Chest color and social status in male geladas. International Journal of Primatology, 30: 791–806.
•  Bouckaert  R.R. (2010) DensiTree: making sense of sets of phylogenetic tree. Bioinformatics, 26: 1372–1373.
•  Delson  E. (1993) Theropithecus fossils from Africa and India and the taxonomy of the genus. In:  Jablonsky  N.G. (ed.), Theropithecus: The Rise and Fall of a Primate Genus. Cambridge University Press, Cambridge, pp. 157–189.
•  Drummond  A.J. and  Bouckaert  R.R. (2015) Bayesian Evolutionary Analysis with Beast. Cambridge University Press, Cambridge.
•  Dunbar  R.I.M. and  Dunbar  E.P. (1975) Contributions to Primatology, Volume 6: Social Dynamics of Gelada Baboons. Karger, Basel.
•  Excoffier  L.,  Laval  G., and  Schneider  S. (2005) Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online, 1: 47–50.
•  Fashing  P.J.,  Nuguyen  N., and  Fashing  N.J. (2010) Behavior of geladas and other endemic wildlife during a desert locust outbreak at Gussa, Ethiopia: ecological and conservation implication. Primates, 51: 193–197.
•  Frost  S.R.,  Jablonski  N.G., and  Haile-Selassie  Y. (2014) Early Pliocene Cercopithecidae from Waranso-Mille (Central Afar, Ethiopia) and the origins of the Theropithecus oswaldi lineage. Journal of Human Evolution, 76: 39–53.
•  Fu  Y.X. (1997) Stastical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics, 147: 915–925.
•  Gamachu  D. (1977) Aspects of Climate and Water Budget in Ethiopia. Addis Ababa University Press, Addis Ababa.
•  Gippoliti  S. (2010) Theropithecus gelada distribution and variations related to taxonomy: history, challenges and implications for conservation. Primates, 51: 291–297.
•  Groves  C.P. (1972) Phylogeny and Classification of Primates. In:  Fiennes  R.N. (ed.), Pathology of Simian Primates. Karger, Basel.
•  Hamilton  A.C. (1982) Environmental History of East Africa: A Study of the Quaternary. Academic Press, London.
•  Hapke  A.,  Zinner  D., and  Zischler  H. (2001). Mitochondrial DNA variation in Eritrean hamadryas baboons (Papio hamadryas hamadryas): life history influences population genetic structure. Behavioral Ecology and Sociobiology, 50: 483–492.
•  Harpending  R.C. (1994) Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Human Biology, 66: 591–600.
•  Heuglin  T., von (1863) Beiträge zur Zoologie Afrikas. Über einige Säugethiere des Bäschlo-Gebietes. Nova Acta Academiae Caesarea Leopoldio, Carolinae Naturae, 30: 1–14.
•  Iwamoto  T.,  Mori  A.,  Kawai  M., and  Bekele  A. (1996) Anti-predator behavior of gelada baboon. Primates, 37: 389–397.
•  Jablonski  N.G. (1993) The phylogeny of Theropithecus. In:  Jablonsky  N.G. (ed.), Theropithecus: The Rise and Fall of a Primate Genus. Cambridge University Press, Cambridge. pp. 209–224.
•  Jolly  C.J. (1972) The classification and natural history of the Theropithecus (Simopithecus) (Andrews, 1916), baboons of the African Plio-Pleistocene. Bulletin of the British Museum (Natural History) Geology, 22: 13–24.
•  Kawai  M. (1979) Contributions to Primatology, Vol. 16: Ecological and Sociological Studies of Gelada Baboons. Karger, Basel.
•  Kifle  Z.,  Belay  G., and  Bekele  A. (2013) Population size, group composition and behavioural ecology of geladas (Theropithecus gelada) and human–geladas conflict in Wonchit Valley, Ethiopia. Pakistan Journal of Biological Science, 16: 1248– 1259.
•  Kimura  M. (1980) A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution, 16: 111–120.
•  Larkin  M.A.,  Blackshields  G.,  Brown  N.P.,  Chenna  R.,  McGettigan  P.A.,  McWilliam  H.,  Valentin  F.,  Wallace  I.M.,  Wilm  A.,  Lopez  R.,  Thompson  J.D.,  Gibson  T.J., and  Higgins  D.G. (2007) Clustal W and Clustal X version 2.0. Bioinformatics, 23: 2947–2948.
•  Matsubayashi  K.,  Shotake  T., and  Belay  G. (1998) Hematological values and plasma erythropoiechin level of wild gelada baboons. Primate Research, 14: 273 (in Japanese).
•  Mittermeier  R.A.,  Rylands  A.B., and  Wilson  D. (2013) Handbook of the Mammals of the World, Vol. 3: Primates. Lynx Edicions, Barcelona.
•  Mori  A. and  Belay  G. (1990) The distribution of baboon species and a new population of gelada baboons along the Wabi- Shebeli river, Ethiopia. Primates, 31: 495–508.
•  Mori  A.,  Iwamoto  T., and  Bekele  A. (1997) A case of infanticide in a recently found gelada population in Arsi, Ethiopia. Primates, 38: 79–88.
•  Mori  A.,  Iwamoto  T.,  Mori  U., and  Bekele  A. (1999). Sociological and demographic characteristics of a recently found Arsi gelada population in Ethiopia. Primates, 40: 365–381.
•  Mori  A.,  Belay  G., and  Iwamoto  T. (2003) Changes in unit structure and infanticide observed in Arsi geladas. Primates, 44: 217–223.
•  Nei  M. (1972) Genetic distance between populations. American Naturalist, 106: 283–292.
•  Nei  M. (1987) Molecular Evolutionary Genetics. Columbia University Press, New York.
•  Rogers  A.R. and  Harpending  H. (1992) Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology and Evolution, 9: 552–569.
•  Rüppell  W.P.E.S. (1835) Neue Wirbelthiere zu der Fauna von Abbyssinien gehörig. Säugethiere. Siegmund Schmerber, Frankfrut-am-Main.
•  Saitou  N. and  Nei  M. (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution, 4: 406–425.
•  Schneider  S. and  Excoffier  L. (1999) Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rate vary among sites: application to human mitochondria DNA. Genetics, 152: 1079–1089.
•  Shotake  T. (1981) Population genetical study of natural hybridization between Papio anubis and P. hamadryas. Primates, 22: 285–308.
•  Shotake  T. and  Nozawa  K. (1984) Blood protein variations in baboons II. Genetic variability within and among herds of gelada baboon in the central Ethiopian plateau. Journal of Human Evolution, 13: 265–274.
•  Snyder-Mackler  N.,  Alberts  S.C., and  Bergman  T.J. (2014) The socio-genetics of a complex society: female gelada relatedness patterns mirror association patterns in a multilevel society. Molecular Ecology, 23: 6179–6191.
•  Tajima  F. (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics, 123: 585–595.
•  Tamura  K.,  Peterson  D.,  Peterson  N.,  Stecher  G.,  Nei  M., and  Kumar  S. (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolution, 28: 2731–2739.
•  Wolde  M. (1970) An Atlas of Ethiopia. Il Poligrafico Priv. Ltd. Co., Asmara.
•  Yalden  D.W. and  Largen  M.J. (1992) The endemic mammals of Ethiopia. Mammal Review, 22: 115–150.
•  Yalden  D.W.,  Largen  M.J., and  Kock  D. (1977) Catalogue of the mammals of Ethiopia. Primates. Monitore Zoologico Italiano (Nuova Serie) Supplemento, 9: 1–52.
•  Yalden  D.W.,  Largen  M.J.,  Kock  D., and  Hillman  J.C. (1996) Catalogue of the mammals of Ethiopia and Eritrea. 7. Revised checklist, zoogeography and conservation. Tropical Zoology, 9: 73–164.
•  Zachos  F.E.,  Apollonio  M.,  Bärmann  E.V.,  Festa-Bianchet  M.,  Gohlich  U.,  Habel  J.C.,  Haring  E.,  Kruckenhauser  L.,  Lovari  S.,  McDevitt  A.D.,  Pertoldi  C.,  Rössner  G.E.,  Sánchez-Villagra  M.R.,  Scandura  M., and  Suchentrunk  F. (2013) Species inflation and taxonomic artefacts—A critical comment on recent trends in mammalian classification. Mammalian Biology, 78: 1–6.