Phylogenetics, genome diversity and origin of modern leopard, Panthera pardus

Leopards, Panthera pardus, are widely distributed across southern Asia and sub‐Saharan Africa. The extent and phylogeographic patterns of molecular genetic diversity were addressed in a survey of 77 leopards from known geographical locales representing 13 of the 27 classical trinomial subspecies. Phylogenetic analysis of mitochondrial DNA sequences (727 bp of NADH5 and control region) and 25 polymorphic microsatellite loci revealed abundant diversity that could be partitioned into a minimum of nine discrete populations, tentatively named here as revised subspecies: P. parduspardus, P. p. nimr, P. p. saxicolor, P. p. fusca, P. p. kotiya, P. p. delacouri, P. p. japonensis, P. p. orientalis and P. p. melas. However, because of limited sampling of African populations, this may be an underestimate of modern phylogeographic population structure. Combined phylogeographic and population diversity estimates support an origin for modern leopard lineages 470 000–825 000 years ago in Africa followed by their migration into and across Asia more recently (170 000–300 000 years ago). Recent demographic reductions likely have led to genetic impoverishment in P. p. orientalis and in the island subspecies P. p. kotiya.


Introduction
The leopard, Panthera pardus , one of the most widely distributed and adaptable big cats, has pelage hues that vary from pale yellow to deep golden or tawny, and are patterned with black rosettes. Melanistic forms occur throughout its range, mostly in humid areas (Seidensticker & Lumpkin 1991;Nowell & Jackson 1996). The coat and colour patterns vary widely across various types of habitat. Pocock (1932) described four different colouration patterns that correspond to semidesert, savannah, rain forest and high mountain leopards. In the Russian Far East the leopard inhabits snowy temperate forests with winter temperatures reaching -25 ° C, and displays a pale cream-coloured long-hair winter coat that has led to its confusion with the snow leopard, Panthera uncia (Pocock 1930). Leopards occur at sea level (Africa, Arabia, India, Java), in foothill areas, in mountains, and on tops of volcanoes (Morocco, Turkmenistan, Iran, Russia, Java). The leopards are found in the Himalayas where they are sympatric with snow leopards up to 5200 m. Throughout their range the leopard feeds on a broad range of prey, including small rodents, birds, different species of ungulates and livestock (Hoogerwerf 1970;Nowell & Jackson 1996;Christen 2000).
The leopard's historic range spanned all of the sub-Saharan and north Africa, the Middle East and Asia Minor, South and Southeast Asia, and extended to the Amur Valley in the Russian Far East. Island ranges included Sri Lanka, Java, Zanzibar, and Kangean (Seidensticker & Lumpkin 1991;Nowell & Jackson 1996). Leopards are still found throughout most of their historic range (Fig. 1), although their numbers have been significantly reduced over the last hundred years due to increasing human population expansion, habitat loss, hunting, and poaching. In some Map showing: (i) historic (dark and light grey) and present (dark grey) geographical distribution of leopards; (ii) distribution of named classical leopard subspecies (big and small three-letter codes together) and distribution of revised subspecies classifications based on present work extending that of Miththapala et al. (1996) (big three-letter codes); and (iii) sample collection sites and number of samples from each site (circles). Full Latin trinomial subspecies names, subspecies codes, countries of their origins and collaborators that provided leopard samples are shown in Table 1. (1) P. p. pardus in Africa; (2) P. p. saxicolor in central Asia; (3) P. p. fusca in India; (4) P. p. kotiya in Sri Lanka; (5) P. p. melas in Java; (6) P. p. orientalis in Russian Far East; (7) P. p. japonensis in North China; and (8) P. p. delacouri in South China. The distinctiveness of the East Asian subspecies was not well supported due to limited sampling, and the molecular taxonomic definition of these subspecies remains provisional (Miththapala et al . 1996).
In the present study we revisit the assessment of molecular genetic variation and genetic differentiation in contemporary leopard populations. We examined leopard geographical partitioning with 25 feline-specific microsatellites, or short tandem repeat loci (STRs), and DNA sequence variation in two mtDNA regions, part of the NADH-5 gene (611 bp) and the control region (CR, 116 bp). Our samples included 36 specimens from Miththapala et al . (1996) and an additional 41 specimens from other geographical locations (Table 1, Fig. 1). We used new methods and additional samples to test for subspecies/population differentiation and to compare the amount of genetic variation within identified leopard subspecies. The results were interpreted in terms of evolutionary history and phylogeography of the leopard in its natural habitat.

Samples
Seventy-seven samples from leopards of known geographical origin were used (Table 1, Fig. 1). Subspecies Panthera pardus saxicolor , P. p. delacouri , P. p. japonensis , and P. p. melas were represented by captive bred individuals only. Samples were selected from unrelated (to the best of our knowledge) leopards with the exception of two P. p. melas which appeared to be relatives (A. Shoemaker, personal communication). Both traditional and revised subspecies names for each sampled leopard are listed in Table 1. Twentytwo wild tiger samples from five subspecies were used as the outgroup for microsatellite analysis (Wentzel et al . 1999). Individual samples of tigers, lions, jaguars, and snow leopards were used as outgroup species in mtDNA analysis. DNA was extracted from whole blood, white blood cells, or fibroblast cultures from skin biopsies using a standard phenol-chloroform method (Sambrook et al . 1989). DNA from plasma was extracted using QIAamp DNA Blood Midi Kits (QIAGEN).
Polymerase chain reactions (PCR) (25 µ L) were performed using 2.5 ng of genomic DNA in 10 m m Tris-HCl (pH 8.3), 50 m m KCl, 1.5 m m MgCl 2 , 200 µ m each of dATP, dCTP, cTTP, dGTP, 1 µ m of each primer and 1 unit Taq-Gold DNA polymerase. For each reaction 35 cycles were performed with 0.5 min denaturation at 94 ° C, 1.5 min annealing at 50 ° C for NADH-5 gene and 55 ° C for CR, and 1 min extension at 72 ° C. Products were checked in 2% agarose gel in TBE buffer. PCR products were purified with CENTRICON-100 filters (Amicon).
The NADH-5 segment was sequenced in both forward and reverse directions using an ABI BigDye Terminator kit. The CR fragment was sequenced twice in the forward direction using an ABI FS Dye Primer kit (in this case PCR primers were designed to include M13 tails). Products were sequenced using Applied Biosystem 373 A and 377 Automated DNA sequencers. Sequences are deposited in GenBank (accession numbers AY035227-AY035292).
Sixty-nine from the total 77 leopards were taken to final analysis (only those successfully amplified and sequenced for total length of both mtDNA segments were taken). Preliminary results from separate analyses of NADH-5 sequences produced phylogenetic associations similar to those when two segments were combined together [shown in Fig. 5 in form of linearized tree derived from neighbourjoining (NJ) tree topology]. The resolution of the maximum parsimony (MP) tree, however, was higher when both mtDNA segments, NADH-5 and CR, were analysed together. The probability of estimating the correct tree seems to be higher when data from different genes are combined in the phylogenetic analysis Huelsenbeck et al . 1996). Thus, we present our final data from combined sequenced analysis. To correct for different mutation rates in NADH-5 and CR, the among site variation option (Yang & Kumar 1996) was applied throughout all mtDNA analysis.
Lion ( P. leo ), tiger ( P. tigris ), jaguar ( P. onca ), and snow leopard ( P. uncia ) individuals were sequenced for the homologous segments of NADH-5 and CR to be used as outgroup species in the phylogenetic analysis. CR was not successfully amplified in lions using primers reported here; thus, only NADH-5 sequences were used.
Of 77 leopard samples, 75 were included in the microsatellite analysis. A DNA sample of a P. p. nimr (Ppa 89) extracted from a museum pelt that amplified only for 10 of the loci and DNA from a Javan leopard plasma sample (Ppa 195) with low yield were excluded from the analysis.
PCR amplifications for each microsatellite locus were performed as described in Menotti- Raymond et al. (1999). PCR amplification products were diluted with sterile deionized water in individual tubes, then multiplex pooled into groups of 4 -5 loci based on product size and fluorescent dye label. Products were resolved by electrophoresis in an ABI Prism TM 310 Genetic Analyser. Data were analysed using ABI Prism TM genescan 2.1 and genotyper 2.0. Amplification products of at least two samples were electrophoresed in each run as standards to correct allele size differences if necessary. PCR product length was used as an actual repeat length (Ellegren 1995).

Data analysis
mtDNA. Leopard NADH-5 and CR sequences were edited and aligned using sequencher (Gene Codes Co., Ann Arbor, Michigan) and clustalx (Thompson et al. 1997) and checked visually. The NADH-5 mtDNA gene sequences were translated into 203 codons. Thirteen sites in leopard sequences could not be unambiguously scored and were excluded from the analysis.
Phylogenetic relationships among mtDNA haplotypes were assessed using three methods implemented in paup* version 4.0 (Swofford 1998): maximum parsimony (MP), minimum evolution (ME), and maximum likelihood (ML). The MP analysis was conducted using a heuristic search, with random addition of taxa and tree-bisectionreconnection (TBR) branch swapping. In the ME approach, neighbour-joining (NJ) trees (Saitou & Nei 1987) were generated with Kimura 2-parameter (Kimura 1980) γcorrected distances; TBR branch swapping was used next to find a minimum evolution tree. The shape parameter (α) for the γ-distribution was estimated using the pamp program in the paml software package (Yang & Kumar 1996;Yang 1999), and was equal to 0.29. The ML analysis was performed using the HKY85 model (Hasegawa et al. 1985) with the among site variation option (α set to 0.29 as estimated by paml). Reliability of all trees was tested by 100 bootstrap replications (Hillis & Bull 1993).
Five different scenarios of leopard geographical subdivision were tested by F ST (with γ-corrected Kimura 2parameter distances) using the amova algorithm (Excoffier et al. 1992) as implemented in the arlequin 1.1 (Schneider et al. 1997). P. p. nimr and P. p. melas were excluded from subdivision analysis due to limited sample size. Parameters of genetic variability for leopard populations were assessed with mega 1.01 (Kumar et al. 1993) and arlequin 1.1, and were measured in terms of polymorphic sites, number of pairwise differences, and nucleotide diversity (genetic diversity of populations with a sample size of four or more was estimated).
The approximate age of modern leopard lineages was estimated using lintre (Takezaki et al. 1995). This program tests the molecular clock on a given topology of a phylogenetic tree and makes linearized trees re-estimating branch lengths under the assumption of a constant rate of evolution (Takezaki et al. 1995). A phylogenetic tree was constructed using the NJ method (Saitou & Nei 1987) and Kimura 2-parameter γ-corrected distances for NADH-5 sequences (611 bp) only. The parameter α for γ-correction for the NADH-5 sequences was estimated using paml (α = 0.90). Both the two-cluster and branch-length tests implemented in lintre were applied (Takezaki et al. 1995). The coalescence point between leopard and lion haplotypes was chosen to be a calibration point and two fossil dates were used. First, 3.5 Ma was used because it has been recorded as the earliest fossils in Tanzania, Africa, for both leopards and lions (Turner & Anton 1997). Second, 2 Ma was used because it is believed to be the lowest bound for the proposed time of split in the Panthera lineage (O'Brien et al. 1987;Wayne et al. 1993). Recent evidence suggests that the snow leopard is a basal divergence in the Panthera genus ( Johnson & O'Brien 1997), therefore the snow leopard was used as an outgroup.
Microsatellites. Pairwise genetic distances among individual leopards using microsatellite data were estimated based on the proportion of shared alleles (Dps) and the kinship coefficient (Dkf ) (Bowcock et al. 1994) with [1 − ps/kf] option as implemented in microsat (Minch et al. 1995). The program neighbor (included in phylip package; Felsenstein 1985a) was used to construct NJ phylogenetic trees from the distance matrices. Nei's D a genetic distances (Nei et al. 1983) computed with dispan (Ota 1993) were used for pairwise comparisons among leopard populations. D a genetic distances are not proportional to evolutionary time, but have been shown to generate correct phylogenetic trees under various evolutionary conditions . One hundred bootstrap iterations were used to estimate the reliability of nodes uniting leopard individuals; one thousand iterations were used for the leopard population trees (Felsenstein 1985b).
Five different subdivision scenarios among leopard populations based on STR data were assessed by R ST , sum of squared size differences (Slatkin 1995), in arlequin 1.1. A Mantel correlation test between pairwise F ST and R ST values was applied with 1000 iterations using mantel 2 (Liedloff 1999) to test whether or not 'subspecies' subdivisions estimated with mtDNA and microsatellite data are congruent. Tests for significance of deviation from Hardy-Weinberg equilibrium for each locus in each population (Guo & Thompson 1992), and tests for genotypic linkage disequilibrium for each pair of loci in each population were performed using genepop software (version 3.1) (Raymond & Rosset 1995). Significant deviations from Hardy-Weinberg equilibrium (α = 0.05) showing deficiency of heterozygotes were found at three loci in three populations: P. p. japonensis (FCA126), P. p. pardus (FCA097), and P. p. kotiya (FCA441) after a Bonferroni correction (Rice 1989).
Variability across 25 microsatellite loci for each leopard population (for those with a sample size of four or more) was measured in terms of percentage of polymorphic loci, average expected heterozygosity, average number of alleles per locus, percentage of unique alleles, average and maximum range of microsatellite repeats, and average variance of microsatellite repeats. Measures of genetic variation were estimated using microsat and exel. Estimates of leopard microsatellite diversity may be biased since only polymorphic loci were used in the present study.

Phylogenetic analysis
To test for evidence of leopard population/subspecies differentiation, the relationships among individual leopards were examined by phylogenetic analyses of mtDNA haplotypes and of composite genotypes from 25 microsatellite loci.
mtDNA. The two combined mtDNA regions of 69 leopards revealed 50 variable sites (44 in NADH-5 and six in the CR), which defined 33 haplotypes in leopards ( Table 2). The haplotypes of each leopard are listed in Table 1. All 13 classical leopard subspecies sampled had private mtDNA haplotypes, i.e. a haplotype found in only one subspecies ( Table 2). The phylogenetic analysis of haplotypes using MP, ME and ML produced concordant topologies (Fig. 2). African and Asian leopards assorted into separate monophyletic groups with two exceptions. The exceptions involved two haplotypes, one found in two Panthera pardus melas (Mel1) leopards and the other in P. p. nimr (Nim1) leopard. These differed from other haplotypes and did not consistently cluster with African or Asian clusters (Fig. 2).
Two clusters of African leopard haplotypes (labelled as PAR-I and PAR-II in Fig. 2) were resolved with relatively high bootstrap support (76% and 84% for the group I and II, respectively, in MP, 78% and 79% in ME, and 87% and 84% in ML). These two groups along with the P. p. nimr (Nim1) were basal in the MP analysis relative to a cluster of Asian leopards (Fig. 2a). The ME tree topology differed in that two clusters of African leopards, PAR-I and PAR-II grouped together (63%) and defined sister taxa with the P. p. nimr haplotype (67%) (Fig. 2b).
Haplotypes from the African PAR-I group were represented by several P. p. shortridgei leopards from different countries of southern Africa (Table 1), and also by single representatives of other African classical subspecies: P. p. panthera (Pan1), P. p. reichenowi (Rei1) and P. p. suahelicus (Sua1) ( Table 1, Fig. 1). The four African individuals with PAR-II haplotypes were restricted to P. p. shortridgei: three from Kruger National Park, South Africa (haplotypes Sho6-Sho8) and one from Zimbabwe (Sho9 ; Table 1).
Microsatellites. Using the composite genotypes of 25 microsatellite loci from 75 leopard individuals, NJ phylogenetic trees were constructed using Dps and Dkf genetic distances with [1 − (ps/kf ) ] option. Both distances produced a similar topology: in all trees leopard individuals tended to cluster together according to their geographical origins, forming eight groups (Fig. 3). In contrast to mtDNA data, there was no evidence from microsatellite genotypes for the PAR-I/PAR-II subdivision among the African leopards. Leopards from mitochondrial groups PAR-I and PAR-II clustered together and were not significantly distinctive with microsatellite R ST value (R ST = 0.009; P = 0.297). The single P. p. sindica associated with P. p. saxicolor leopards, forming a group of central Asian leopards. Sri Lankan leopards (P. p. kotiya) were grouped in the same cluster, but closely aligned with Indian leopards, P. p. fusca (Fig. 3). Among East Asian leopards three classical subspecies © 2001 Blackwell Science Ltd, Molecular Ecology, 10, 2617-2633 Table 2 Haplotypes and variable sites in combined analysis of the NADH-5 (611 bp) and CR (116 bp) mtDNA in leopard (P. pardus) and outgroup species, lions (P. leo), tigers (P. tigris), jaguars (P. onca) and snow leopard (P. uncia) NADH5 CR

Population subdivision and subspecies recognition
To evaluate the extent of population differentiation in leopards we tested five different geographical scenarios and compared them based on analysis of molecular variance (amova) with both mtDNA and microsatellite data. The first scenario considered only two groups: all African leopards (P. p. shortridgei, P. p. panthera, P. p. reichenowi, P. p. suahelicus) vs. all Asian leopards (P. p. saxicolor, P. p. sindica, P. p. fusca, P. p. kotiya, P. p. delacouri, P. p. japonensis, and P. p. orientalis); these major groups were proposed based on the topology of the mtDNA phylogenetic trees (Fig. 2). amova performed using mtDNA found 68.9% of the variation (F ST , Table 3) between the continents and 31.1% within the continents (P < 0.0001). With microsatellites,  (Table 2). Individual samples of Panthera leo, P. tigris, P. onca, and P. uncia are taken as outgroup species. (a) Maximum parsimony (MP) tree. MP tree constructed with paup* (Swofford 1998) and a general heuristic search; with random taxon addition and treebisection reconnection branch swapping; shown in majority rule consensus of 10 145 trees (length = 212, CI = 0.684). Numbers above branches represent bootstrap support (100 replicates); only those with > 50% are shown. Numbers below show number of steps/number of homoplasies. (b) Minimum evolution (ME) tree, constructed with paup* using Kimura 2-parameter γ-corrected distances (α = 0.29 as defined by paml) and neighbour-joining algorithm followed by tree-bisection-reconnection branch swapping. Values above branches represent support from 100 bootstrap replicates (only those with > 50% are shown). Maximum Likelihood (ML) approach performed using HKY85 model (Hasegawa et al. 1985) with the among site variation option (the α set to 0.29 and estimated with palm) produced generally the same topology as ME tree (not shown). Bootstrap support from 100 replicates for ML tree is shown in ME tree (b) below branches (only those with > 50% are shown). Table 3) could be explained by continent subdivision, and 63.9% of the variation was retained within the continents (P < 0.0001). In a second test, leopards were considered as three geographical groups: (1) African (as above); (2) Central Asian, Indian and Sri Lankan together (P. p. saxicolor, P. p. sindica, P. p. fusca, and P. p. kotiya); and (3) East Asian (P. p. delacouri, P. p. japonensis, and P. p. orientalis). This scheme was suggested by the topology of the microsatellite trees (Fig. 3). With this grouping, 63% of mtDNA and 31.6% of microsatellite variation can be explained by geographical partitioning (Table 3), and 37% and 68.4% of the variation, respectively, are retained within the groups (P < 0.0001). The third scenario involved four geographical groups and differed from the second in that central Asian (P. p. saxicolor, P. p. sindica), and Indian and Sri Lankan (P. p. fusca, and P. p. kotiya) were analysed separately. Among group variation in this case was slightly higher than in previous scenario, 63.3% with mtDNA (P < 0.0001) and 32.9% with micro-satellites (P < 0.0001; Table 3).

36.1% of the variation (R ST ,
In the fourth scenario applied, leopards were divided into seven different groups that correspond to the revised subspecies, based on the mtDNA and microsatellite phylogenetic analysis (Figs 2, 3 and 4). With this grouping, 76.04% of mtDNA variation was distributed among leopard populations and only 23.96% (P < 0.0001) within the populations (Table 3). With microsatellite data, 35.8% of the variation was found between the subspecies, while  (Table 1). Trees constructed based on proportion of shared alleles (Dps) and kinship coefficient (Dkf ) genetic distances with 1 − (kf/ps) option in microsat (Minch et al. 1995) produced the identical topologies; Dps tree is shown. Numbers are individual Ppa identification (Table 1).
64.2% of the variation occurred within the subspecies (P < 0.0001). When African leopards were considered as two groups, PAR-I and PAR-II (the last scenario), subdivision was even higher: F ST = 83.8%, and R ST = 36.3 (P < 0.0001; Table 3). Based on these tests, variation among leopard populations was best explained by grouping leopards on the basis of 'subspecies' scenarios (see below). The statistical significance of pairwise population differentiation was tested by F ST for mtDNA data and by R ST for microsatellite data (Table 4). Each population was significantly different from the others by pairwise F ST for mtDNA data. For microsatellites, two population comparisons of R ST , PAR-I and PAR-II, and DEL and JAP, were not significantly different (P = 0.297 and 0.069, respectively). R ST pairwise comparisons performed with African leopards considered to be a single group (PAR) revealed all revised subspecies to be distinctive (P < 0.001). Mantel correlation analysis between F ST and R ST pairwise values (given in Table 4) revealed the two matrices to be significantly correlated (P < 0.005; g = 5.661, Z = 20138 and r = 0.841).
In summary, the position of each mtDNA haplotype (Fig. 2) and each leopard's composite microsatellite genotype in phylogenetic trees (Fig. 3) correlated well with their geographical origins; however, there was not strong bootstrap support for the phylogeographic clusters. The phylogeographic concordance of both mtDNA and microsatellite analyses plus the significant partitions of distinctive groups (Tables 3 and 4) would support the recognition and genetic distinctions for a minimum of nine groups: P. p. pardus- PAR, P. p. saxicolor-SAX, P. p. nimr-NIM, P. p. fusca-FUS, P. p. kotiya-KOT, P. p. delacouri-DEL, P. p. japonensis-JAP, P. p. orientalis-ORI, and P. p. melas-MEL, which we propose to recognize as revised subspecies of P. pardus.
A phylogenetic analysis of seven revised leopard subspecies (i.e. those with multiple individuals) was constructed using D a genetic distances for microsatellite population data (Fig. 4). The deepest split separated African leopards, P. p. pardus (combined PAR-I and PAR-II) from other leopard groups. Central Asian leopards, P. p. saxicolor, followed Africans in the phylogenetic tree. Indian leopards (P. p. fusca) clustered with Sri Lankan (P. p. kotiya) individuals with relatively high bootstrap support (83%). The three East Asian subspecies formed a monophyletic lineage with high statistical support (99%): P. p. orientalis consistently grouped with P. p. japonensis; and the pair formed a sister taxon with P. p. delacouri (Fig. 4). This phylogenetic tree corresponds rather well with the geographical distribution of leopard populations. ¶The same as 7 groups, but African divided into P. p. pardus I (PAR-I) and P. p. pardus II (PAR-II). P. p. nimr and P. p. melas were excluded from all scenarios due to limited sampling. †Calculated with Kimura 2-Parameter distances (Kimura 1980). ‡Calculated with R ST option in arlequin 1.1 (Schneider et al. 1997).  (Nei et al. 1983) genetic distances. Numbers indicate bootstrap percentages.

Genetic variation
Estimates of population genetic variability calculated from mtDNA sequences and microsatellite loci for each revised leopard subspecies are summarized in Table 5. Genetic diversity varied appreciably among the leopard subspecies, and trends were similar for both mtDNA sequences and microsatellites. Estimates of genetic variation were the highest in African leopards, P. p. pardus, lowest in the Far Eastern leopard P. p. orientalis, substantially reduced in Sri Lankan P. p. kotiya leopards, and moderate in other populations. With mtDNA data, for example, there was 21 variable sites in P. p. pardus population while there was only one in P. p. orientalis or P. p. japonensis. Mean number of pairwise nucleotide differences in P. p. pardus (8.77) was more than 50 times higher than within P. p. orientalis (0.17) or about 15 times higher than in the P. p. kotiya (0.56). Nucleotide diversity (π) was high in P. p. pardus (1.22%), very low in P. p. orientalis (0.02%), and in between these two extremes in other subspecies populations (Table 5). The African leopards, P. p. pardus, also had the highest diversity estimates for microsatellite loci. Heterozygosity in this population was the highest (0.803), somewhat lower in P. p. fusca (0.696) and in P. p. saxicolor (0.616), relatively low in P. p. kotiya (0.485), and lowest in P. p. orientalis (0.356). The average number of microsatellite alleles, range of microsatellite repeats and microsatellite variance had the same trends: highest in P. p. pardus (8.52, 9.72 and 7.28, respectively), lowest in P. p. orientalis (2.60, 2.84 and 1.71) and moderate in all others (Table 5). The microsatellite allele size distribution was most heterogeneous in African leopards (allele size distributions for each population and each locus are given as complementary information to this paper at http://lgd.nci.nih.gov). Alleles found in P. p. orientalis and in P. p. kotiya populations were almost always a subset of those seen in P. p. japonensis and P. p. fusca, respectively, and they were discontinuously distributed, which may suggest a founder effect in history of P. p. orientalis and P. p. kotiya populations.

Diagnostic characteristics
Each revised leopard subspecies possessed populationspecific mtDNA haplotypes, and/or microsatellite alleles (Table 6). Diagnostic mitochondrial sites were found in every revised subspecies except P. p. japonensis. Two mtDNA sites were specific for P. p. orientalis; one for P. p. delacouri, P. p. fusca and P. p. kotiya; three sites were unique for P. p. melas and five for P. p. nimr (Tables 2 and 6). Three fixed sites present in P. p. pardus leopards are not considered to be diagnostic for the African group since these sites were found in P. p. nimr as well, suggesting close evolutionary relatedness of these two subspecies (Table 2). Three sites were shared among P. p. pardus, P. p. nimr and P. p. melas leopards as well as with most outgroup species (P. tigris, P. leo, P. onca, and P. uncia) ( Table 2). African leopards in general had the largest number of mitochondrial sites in common with outgroup species (Table 2).
All revised subspecies, except P. p. saxicolor and P. p. melas, revealed subspecies-specific microsatellite alleles ( Table 6). Frequencies of such private alleles, however, were low in each population (2.83 -5.80% of total number of alleles), with the exception of African leopards P. p. pardus, where they were 29.0% (Table 5). Number of diagnostic mtDNA sites and percentage of subspecies-specific microsatellite alleles should be considered relative to the sample sizes of leopard populations presented here.

Estimation of divergence times
The mtDNA sequence divergences were used to test the hypothesis of a molecular clock for different leopard haplotypes and to estimate approximate times of leopard divergence (Fig. 5, Takezaki & Nei 1996 Table 5 Genetic variation across mtDNA gene segments 611 bp,and Control Region,116 bp) and 25 microsatellite loci in seven revised leopard subspecies  124,208; Kot2,Del2, test did not reveal significant rate heterogeneity among leopard sequences (confidence probability; CP < 95%); with the branch-length test two sequences (Del 1 and Del 2) have evolved faster than the rest (for both CP = 95%). However, because these two sequences were only marginally deviated in the branch-length test, they were not excluded, and thus all sequences were used to construct a linearized tree (Fig. 5). Two dates were estimated, the divergence time of the first major node (node 2), and the time of origin of Asian leopards (node 3). Using the equation, H = µT, H was equal to the linearized heights estimated by the two-cluster test, and µ was a substitution rate estimated based on known fossil dates (T). When T was set to be equal to 3.5 million years (Myr) as the time when leopards and lions split (node 1), µ was estimated to be 0.

Discussion
A population genetic and phylogeographic assessment of leopards sampled from specific geographical origins throughout their current range was determined for 77 leopards representing 13 of 27 named classical subspecies. Using DNA sequence data from two mitochondrial gene segments, NADH-5 and CR, and composite microsatellite genotypes for 25 microsatellite loci, evidence for nine revised subspecies was obtained ( Fig. 1 and Table 6). Recognition of the nine revised subspecies was derived from phylogenetic analyses of mtDNA haplotypes (ME, MP and ML algorithms, Fig. 2) plus statistically significant distinct F ST measures (Table 4). The mtDNA distinctions were supported by phylogenetic analysis of composite microsatellite genotypes that assorted individuals into the same groups (Fig. 3). Seven revised subspecies for which there were four or more representatives (all except Panthera pardus melas and P. p. nimr) consistently showed genetic differentiation based on R ST estimates for microsatellite allele distributions among the population (Tables 3 and 4). Diagnostic mtDNA sites, haplotypes or subspecies specific microsatellite alleles form the basis for recognition of the nine subspecies (Table 6). The results affirm and extend the provisional subspecies categories proposed earlier (Miththapala et al. 1996), based on other molecular evolutionary markers. The two revised subspecies with the fewest sampled individuals, P. p. melas (n = 2) and P. p. nimr (n = 1), should be considered as tentative at this stage, although highly distinctive mtDNA haplotypes were observed in multiple individuals from these regions. Populations of P. p. nimr appear to have been isolated for quite a long time, accumulating multiple diagnostic sites that distinguish it from any other subspecies (Table 2). Presently, not more than 200 leopards are thought to be left in the whole Arabian peninsula (Lagrot & Lagrot 1999).
Sampling in central, west and northern Africa was limited (one P. p. suahelicus, one P. p. reichenowi, and one P. p. panthera, Fig. 1 and Table 1). These isolated individuals did not discriminate from the more numerous southern African samples (Figs 2 and 3). More extensive sampling in the future may reveal further partitions among north/ central African leopards. The occurrence of divergent mtDNA haplotype lineages (PARI and PARII, Fig. 2) may reflect ancestral subdivisions consistent with this possibility. Nonetheless, until better evidence is developed, we consider African leopards to comprise a single revised subspecies, P. p. pardus.
In east Asia, we found significant differentiation and diagnostic markers that support the recognition of P. p. orientalis, P. p. japonensis and P. p. delacouri. Miththapala et al. (1996) had correctly deferred judgement on these three subspecies due to inadequate sampling, which is remedied here. The island population of Sri Lanka, P. p. kotiya, was distinctive, but was closely aligned to the mainland India subspecies P. p. fusca. This similarity is likely due to an historic origination of P. p. kotiya from P. p. fusca founders. The Javan leopard P. p. melas was highly distinctive from other Asian leopards for evolutionary reasons that remain uncertain (see Miththapala et al. 1996).

Genetic diversity
The amount of genetic diversity revealed in the leopards was comparable to or higher than those reported for other cat species, such as lions, cheetahs (Driscoll 1998), jaguars (Eizirik et al. 2001), andpumas (Culver et al. 2000). The genetic variation in leopards, however, varied significantly across their geographical range. There may be bias in estimation of actual genetic diversity among leopard populations due to different and sometimes insufficient sample sizes. Three populations (P. p. saxicolor, P. p. japonensis and P. p. delacouri) were represented only by captive-bred individuals, and thus, the reported genetic estimates reflect the status of captive populations of these revised subspecies.
The African leopards were the most genetically variable among all leopard subspecies by both mtDNA sequences and microsatellites (Table 5). Similar high diversity has been reported for Tanzanian leopards (Spong et al. 2000). The lowest level of genetic variation in both types of markers was observed in the Far Eastern leopard, P. p. orientalis. This population has a documented history of demographic and range reduction, and it is the most critically endangered leopard subspecies (Miquelle et al. 1996;Nowell & Jackson 1996;Uphyrkina et al., unpublished data). Sri Lankan leopards, P. p. kotiya, which had previously been reported as showing diminished genetic diversity with several genetic metrics (Miththapala et al. 1991(Miththapala et al. , 1996, also showed relatively low levels of microsatellite variation (Table 5). The mtDNA diversity was somewhat low in P. p. saxicolor (mean number of pairwise differences was 0.50 and nucleotide diversity was 0.07%). This may be explained by the sampling of only captive individuals which are thought to be substantially inbred (A. Shoemaker, personal communication).

Radiation of modern leopards
African leopards were the first group to split off from the phylogenetic tree based on the microsatellite data (Fig. 4) and they were also basal relative to Asian leopards in the MP phylogenetic tree with mtDNA data (Fig. 2). African leopards possessed the broadest range of genetic variation by all molecular genetic techniques applied to date: allozymes, mtRFLP, minisatellites (Miththapala et al. 1996), microsatellites and mtDNA (Table 5). Further, P. p. pardus share more mitochondrial sites in common with outgroup species than other subspecies ( Table 2). We interpret these observations as indicative of an African origin for leopard genetic diversity which we estimate as between 470 000 and 825 000 years ago depending on which fossil calibration dates were employed (Fig. 5).
The Asian lineages are estimated as somewhat younger, between 170 000 and 300 000 years ago, consistent with a migration out-of-Africa to the middle east and east to eastern Asia during that interval. The leopard may have had to cross the Afro-Arabian landbridge, perhaps by the Egyptian-Sinai-Israeli passageway, following the invasion of many sub-Saharan animal and plant forms into the eastern Mediterranean and Eurasia during late Pliocene and early Pleistocene (Tchernov 1988). The leopard's migration would correspond precisely in time with the postulated migration of modern human populations out-of-Africa similarly deduced from patterning of mitochondrial and nuclear genomic diversity (Hedges et al. 1992;Nei & Roychoudhury 1993;Bowcock et al. 1994;Goldstein et al. 1995;Calafell et al. 1998;Ingman et al. 2000).
Based on fossil dating of the earliest members of Panthera, Hemmer (1976) proposed that first the jaguarlike ancestor of this subgenus spread over Africa, Europe, southern and northern Asia and North America in the middle Lower Pleistocene. Then differentiation into the living species took place as a second stage somewhere in Lower Pleistocene. The ancestral leopard may have gone extinct during faunal turnovers throughout the world except in Africa, and the modern leopard may then have spread out of Africa again. The leopard appears to have taken the same routes that were used by modern human migrations (Hedges 2000).

Conclusions
In the present paper we have confirmed and extended the phylogenetic discrimination of seven phylogeographic groups of leopards to nine revised subspecies, one African, Panthera pardus pardus, and eight Asian subspecies, P. p. saxicolor, P. p. fusca, P. p. kotiya, P. p. melas, P. p. delacouri,