Unravelling the maternal evolutionary history of the African leopard (Panthera pardus pardus)

Data set preparation

This study leverages mitochondrial NADH-5 data from all previous African leopard genetic studies to shed light on the evolutionary history of the leopard at the southernmost end of its range. Although the NADH-5 locus has been sequenced in 139 leopard specimens collected in South Africa, these data have never been analysed in combination with all the available data from the rest of Africa. Yet, such an analysis would be required to resolve the geographic distribution, and regions of potential overlap, of the ancient leopard clades PAR-I and PAR-II. This is particularly important, since the Morris et al. (2023) data set also represents the only leopard samples from the massive central plateau, known locally as the Highveld, and from along the entire eastward escarpment transition zone down to the Lowveld. This sample of 33 specimens collected in Mpumalanga, therefore, represents a number of ecologically diverse habitats in a region where two divergent maternal clades/nuclear populations have previously been found to overlap (Morris et al., 2023).

The genetic data set here assembled comprises the largest geographical sample for any African leopard study to date. These include 133 individual leopard specimens from Ropiquet et al. (2015) (South Africa and Mozambique), 56 from Anco et al. (2018) (samples sourced from all over Africa, comprised of a combination of museum specimens and samples collected from live animals by other researchers), 13 from Uphyrkina et al. (2001) (southern Africa–KNP, Botswana, Namibia), nine from Paijmans et al. (2018) (data from historical leopard museum specimens kept at the Natural History Museum of Denmark) and data for 33 specimens collected in Mpumalanga, sampled across the highveld-lowveld gradient of north-eastern South Africa (Morris et al., 2023) (Table 1). The accession numbers for the 33 NADH-5 leopard samples generated from the Mpumalanga province are OQ132962–OQ132992. Almost all studies provided an exact sample collection location, however the museum samples in Paijmans et al. (2018) and Uphyrkina et al. (2001) only specified a country and/or general region within a country of origin. We used a total of 244 NADH-5 sequences in this study, representing leopards from 18 countries across the African continent (Table 1). Consensus sequences for each individual were aligned in BioEdit (v7.0.5.3) (Hall, 1999) and checked for any errors, large gaps or missing data. The 244 samples in the final data set were trimmed to a length of 425 bp, corresponding to positions 12,632–13,058 in the complete mtDNA genome of Felis catus (Lopez, Cevario & Obrien, 1996) (Accession number U20753).

Genetic structure

To explore how mitochondrial DNA variation was structured and to date distinct clades of African leopard, phylogenetic trees were reconstructed from NADH-5 data using the program BEAST (Bayesian Evolutionary Analysis by Sampling Trees (v2.5.1)) (Bouckaert et al., 2019). A Bayesian framework allows the definition of parameters such as mutation rates, clock rates and ancestral calibration points as distributions and not simply point estimates. The substitution model (HKY + G) was selected by jModelTest (v2.1.10) (Darriba et al., 2012) using Bayesian Information Criterion (BIC). Heterogeneity in mutation rates was modelled using a gamma distribution with four bins. The clock model was defined by first running the analysis with relaxed lognormal and exponential clock options. The posterior distribution of the standard deviation of the clock rate did not include zero for the relaxed lognormal clock, suggesting the data fitted this model better than a strict clock. A relaxed lognormal clock rate of 0.0122 was set based on the known NADH-5 mutation rate in other feline species (Culver et al., 2000; Mukherjee et al., 2010).

Two nodes were calibrated with ancestral Panthera lineages. This necessitated the inclusion of the following outgroup taxa: the lion (P. leo, three sequences: accession numbers KP001498, KP001502 and KP001506) and three non-African leopard subspecies (P. p. nimr, P. p. fuscia and P. p. orientalis: Accession numbers AY035279, EF199743 and KY866876). Calibration priors were set up as follows: (1) the divergence between lion and leopard occurring at a mean of 2.57 Ma, with a normal distribution from 1.89 to 3.28 Ma (Figueiro et al., 2017) and (2) the divergence of Asian and African leopards at approximately 710 Ka with a normal distribution from 457–956 Ka. The analysis was run ten times independently with each run consisting of 200,000,000 Markov Chain Monte Carlo (MCMC) generations, sampling trees and parameters every 100,000 generations. Each run log was checked in Tracer (Rambaut et al., 2018). A burn in of 20% was discarded for each run. All post-burn in trees in all runs were combined using Logcombiner and a consensus tree was calculated using maximum clade credibility in Tree Annotator (both programs are in the BEAST package v2.5.1 (Bouckaert et al., 2019)). The above analysis was conducted separately on the full NADH-5 data set of 244 individual African leopards and then on the 47 leopard haplotypes that we observed in the data. Trees were edited in Figtree (v1.4.3) and Treegraph2 (v2.15.0) (Stover & Muller, 2010). Only nodes supported by a posterior probability of greater than 0.5 were shown in the haplotype tree, whereas the posterior probability was given for each node.

To investigate any non-tree-like evolution in African leopards, a median-joining haplotype network was reconstructed in Popart (v1.7) (Leigh & Bryant, 2015). Unlike a Bayesian analysis that can account for poor quality sequences, network reconstruction can be biased by samples with too much missing data. Therefore, sequences containing more than 5% missing data were excluded (28 samples from our dataset, see Table S1), resulting in the omission of five haplotypes compared to the haplotype tree.

Mitochondrial phylogeography

Although some previous studies have shown that PAR-II occurred mainly in southern Africa (Uphyrkina et al., 2001; Anco et al., 2018), its exclusive occurrence only in the southern and south-eastern parts of South Africa and the increasing occurrence of PAR-I in northern South Africa suggests this country as the likely stronghold for PAR-II. We demonstrated this graphically by plotting haplotype distribution and frequencies onto the geographic space of the leopard’s range in Africa.

We further tested this “South Africa” hypothesis using a hierarchical analysis of molecular variance (AMOVA), analysing the haplotype frequencies of 16 hypothetical two-group partitions. These ranged from a scenario where only South Africa haplotypes comprised the first group, and all other African haplotypes the second group, to a hypothesis where all countries in southern Africa (South Africa, Namibia, Botswana, Zimbabwe and Namibia) comprised the first group. All combinations of southern African countries were tested, with South Africa always featuring in the first group. Variance components of the data and resulting F-statistics were tested for significance against a null hypothesis of zero variance using 10,000 permutations.

Diversity

We estimated the genetic diversity of African leopard mitochondrial DNA clades using Arlequin (v3.5.2.2) (Excoffier & Lischer, 2010). Initially we performed the analysis with all the samples in the study (all of Africa) and then determined diversity indices for groups based on the two major clades determined from the phylogenetic tree (PAR-I and PAR-II lineages). To explore the diversity exclusive to leopards in South Africa, we also estimated diversity indices for the following two geographical groups: (1) South Africa and (2) the rest of Africa.

Inferring historical demography

Despite the availability of many leopard mtDNA sequences, a demographic analysis has never been carried out. First, we tested our alignment against a null hypothesis of constant population size using Tajima (1989), Fu (1997) in Arlequin (v3.5.2.2). Positive values for Tajima’s D and Fu’s Fs can denote a past population contraction whereas negative values may indicate a population expansion. We also performed a mismatch distribution in Arlequin by calculating the sum of squared deviations (SSD) and raggedness. Significantly low SSD and raggedness can imply that a population underwent a sudden expansion (Rogers & Harpending, 1992).

We also modelled historical effective population size changes using a Bayesian Skyline model, which also leverages the demographic signal in DNA sequence alignments (Drummond et al., 2005). Here, we analysed PAR-I and PAR-II separately, to infer their independent demographic histories, although we concede that both clades contain significant structuring within them. The evolution of genetic structure can also lead to increases in historic effective population sizes, which can be mistaken for a population expansion on a Bayesian Skyline Plot (BSP). BSPs were reconstructed using BEAST and Tracer (Drummond et al., 2005; Bouckaert et al., 2019). Parameters were set using the same clock rate and models as the initial large phylogenetic tree however, each clade had a prior set based on the divergence time estimated for the corresponding node in the phylogenetic tree. MCMC was set at 50,000,000 permutations with a log every 10,000 retained. As with the phylogenetic analysis, MCMC convergence was determined in TRACER.

Phylogenetic structure

The phylogenetic trees (Figs. 1 and S1) highlighted an initial divergence of African leopard samples into PAR-I and PAR-II clades (posterior = 1). We dated this split to 0.7051 Ma (95% HPD [0.4477–0.9632] Ma) (Fig. 1). Both lineages then diversified at about the same time (PAR-I = 393–918 Ka, PAR-II = 305–846 Ka) giving rise to the observed diversity within each clade. PAR-I contained individuals from across Africa, however, only eight supported PAR-I clades contained more than three haplotypes, indicating that, apart from a major clade in western-central Africa (Cameroon and Gabon) and western Africa (Senegal), most haplotypes did not structure into clades with a clear geographic association (Fig. 1). PAR-II on the other hand, included 112 samples, with only five not originating in South Africa. These non-South Africa PAR-II individuals were sampled either in southern or central Africa (Zimbabwe, Mozambique, Zambia, DRC and Burundi). The two most basal PAR-II haplotypes were sampled in Zambia and the DRC, with the more derived PAR-II haplotypes in southern Africa, with the exception of one sample from Burundi with a haplotype identical to nine samples from South Africa. PAR-II was the only lineage sequenced in all leopards sampled in the Western Cape (n = 27), Eastern Cape (n = 11) and Kwa-Zulu Natal (n = 43) provinces of southern and south-eastern South Africa. The remainder of the South Africa PAR-II individuals were all from the north-eastern part of that country (Mpumalanga and the adjacent Kruger National Park), and these localities contained either PAR-I (n = 31) or PAR-II (n = 25), with the former more prevalent in the escarpment (Lydenburg, n = 3) and Lowveld (Manyeleti GR, n = 11; Andover NR, n = 3; Kruger National Park, n = 14) and the majority of the latter sampled on the Highveld (Loskop Dam, n = 13) and eight in the Lowveld (Manyeleti GR, n = 3; Kruger National Park, ncc = 9) (Figs. 1 and S1).

Reticulate structure

The median-joining haplotype network (Fig. 2) supported the structure observed in the phylogenetic tree (Figs. 1 and S1), confirming a single origin for the PAR-II lineage and showing that PAR-I and PAR-II were the first two lineages to diverge from the common maternal ancestor of all African leopards. Both tree and network also suggest that PAR-II is geographically structured internally, with unique haplotypes occurring in the Western-Eastern Cape, in KwaZulu-Natal, Mpumalanga and the Kruger National Park. PAR-I was also structured in western Africa, but less so in central and southern Africa with high levels of haplotype sharing among countries. PAR-I leopards in central and southern Africa appear to all belong to a phylogenetically diffuse set of haplotypes, all closely related to each other and with the most common haplotype (n = 55) centrally related to almost all others (Fig. 2).

Spatial structure

To visualise how maternal genetic variation was partitioned across geographic space, we plotted the distribution of haplotypes by location (Fig. 3). At a continental scale (Fig. 3A) PAR-I was distributed across most of the leopard’s African range from Algeria to northern South Africa. PAR-II, on the other hand occurs from the DRC and Zambia in Central Africa, with frequencies increasing in a southern direction. At the regional level of southern Africa (Fig. 3B), PAR-II haplotypes are at highest frequency in South Africa, occurring with PAR-I in the South African Lowveld (Kruger National Park and Manyeleti GR) and exclusively in leopards from the South African Highveld, KwaZulu-Natal and the Cape.

AMOVA analyses returned the highest possible F_CT (between group variation) and F_ST (variation among populations within groups) values for the two-group comparison separating PAR-I and PAR-II (Table 2). All F_ST values were significantly different from zero, although we concede that small sample sizes for some countries may have affected values at this hierarchical level. Furthermore, owing to a single degree of freedom in all other two-group scenarios, no F_CT value attained significance, although this statistic did vary appreciably (0.01956–0.11072, Table 2) reflecting the relative levels of variance attributable to each scenario. When several geographic scenarios partitioning different combinations of southern African countries with South Africa were tested, we found that the two-group scenario comparing South Africa with all other leopards maximised the values of both F-statistics, with the comparison including South Africa and Zimbabwe returning the second highest values. Interestingly, the above three grouping scenarios were all more likely than one of a single African leopard population (Table 2).

Diversity and demography

Genetic diversity and demography indices were calculated for PAR-I, PAR-II and the whole data set. We also calculated the genetic diversity of all leopard samples collected from South Africa as we were interested to know what effect the existence of both PAR-I and PAR-II would have on the genetic diversity in that country (Table 3). A total of 47 haplotypes were observed among 244 African leopard NAHD-5 sequences (Table 3), and of these, 28 belonged to PAR-I and 18 to PAR-II. Haplotype diversity was high for both PAR-I (0.858) and PAR-II (0.895), but nucleotide diversity was higher in PAR-I. South Africa contained 34 haplotypes compared to the rest of the continent which displayed 36 different haplotypes (Table 3) and with levels of haplotype and nucleotide diversity exceeding either PAR-I and PAR-II, but slightly lower than the African total.

When we tested these groups against a hypothesis of constant population size, both the total sample and PAR-I returned negative values for Tajima’s D and significantly negative values for Fu’s Fs, indicating a population expansion. On the other hand, only Fu’s Fs was significantly negative in PAR-II and Tajima’s D was positive. In the South Africa sample set, which contained both PAR-I and PAR-II, only Fu’s Fs was negative (Table 3). Similarly, mismatch distributions showed evidence for a population expansion only in the PAR-I lineage, with significant SSD and raggedness values (Table 3).

We modelled past changes in effective population size using Bayesian skyline plots. Our results indicate that PAR-I underwent a population expansion beginning around 150 thousand years ago (Ka), followed by a more recent decline within the last 20 Ka (Fig. 4A). PAR-II, however, did not reflect a population expansion, but rather a steady contraction within the last 100 Ka to its lowest size today (Fig. 4B).

The distinctiveness of PAR-I and PAR-II

As in previous studies, we observed an initial divergence between the two African matrilineages PAR-I and PAR-II (Figs. 1 and 2). However, in this study we attempted to date the timing of splits using secondarily derived divergence estimates from a whole genome study on big cats (Figueiro et al., 2017). This analysis indicated a mid-Pleistocene split of an ancestral female African leopard lineage into the lineages PAR-I and PAR-II occurred between 447 and 963 Ka (Fig. 1). When we examined the data in geographic detail, it became clear that PAR-I did not extend to the southernmost leopard populations of the Cape and KwaZulu-Natal regions of South Africa, which were exclusively PAR-II. By analysing the Mpumalanga samples of Morris et al. (2023) in an African context, we show that their East Mpumalanga and West Mpumalanga clades were in fact none other than PAR-I and PAR-II, respectively. Intriguingly, Morris et al. (2023) also showed that East and West Mpumalanga were also differentiated at 17 microsatellite loci providing the tantalising prospect that PAR-I and PAR-II are also nuclear DNA populations. Other nuclear microsatellite studies did not highlight differences among African leopard samples (Uphyrkina et al., 2001; McManus et al., 2015; Ropiquet et al., 2015). In the case of Uphyrkina et al. (2001), only four individuals contained PAR-II haplotypes and none were from north-eastern South Africa where Morris et al. (2023) detected their major nuclear DNA discontinuity between Highveld and Lowveld samples. McManus et al. (2015) on the other hand, only studied populations within the exclusive range of PAR-II. Ropiquet et al. (2015) did not sample Highveld leopards, but nevertheless obtained a large South Africa data set, which included the Lowveld region and Mozambique which is around the southern limit of PAR-I. Although they claim that they did not find evidence for a nuclear population split across their sample set, it is clear in retrospect that the factorial correspondence analysis they conducted on their microsatellite data (Fig. 4 in Ropiquet et al. (2015)) revealed a major nuclear genetic discontinuity between samples mainly from Mozambique which formed a distinct cluster (PAR-I) that intercepted only with samples from the South African Lowveld (represented by the Kruger National Park) and not with samples from the Cape and KwaZulu-Natal (PAR-II). Thus, PAR-I and PAR-II are quite likely also nuclear DNA populations. Alternatively, distinct mtDNA patterns might also evolve if female leopards disperse to a lesser extent than males as has been suggested previously (Bailey, 1993; Fattebert et al., 2015). However, higher male dispersal may not translate into sex-biased gene flow as males tend to suffer higher natural mortality than females (Balme & Hunter, 2004; Daly et al., 2005). Population-scale sequencing of nuclear genomes sampled from southern Africa would help confirm or refute these interpretations.

The leopards of South Africa

Our analyses of spatial structure showed that although PAR-II haplotypes can be found at a lower frequency in Central Africa, their frequency increases in south-eastern Africa (Zimbabwe, Mozambique and north-eastern South Africa), until they make up the only mtDNA lineage present at the southernmost extent of the species distribution (Fig. 3). AMOVA indicated that the highest variation between groups and among populations within groups was between South Africa and all other countries. Taken together the geographic distribution of haplotypes (Fig. 3) and the AMOVA results (Table 2) suggest strongly that the mtDNA lineage PAR-II, and possibly its associated nuclear DNA population, are typical of leopards of the South African Highveld and coastal regions from KwaZulu-Natal to the Cape of Good Hope. This part of South Africa is thus the stronghold for PAR-II diversity. The range of PAR-II up the western coast or in the Kalahari and Namib regions of north-western South Africa, southern Namibia and southern Botswana is yet to be determined (black question marks, Fig. 3B) and would require dense sampling in these areas.

Evolutionary history of the African leopard

We used our data set to explore the maternal evolutionary history of leopards on the African continent. The divergence of PAR-I and PAR-II from a common ancestor occurred during the mid-Pleistocene transition (1,250–700 Ka) (Clark et al., 2006; Legrain, Parrenin & Capron, 2023), when glacial cycles increased from 40 to 100 ka, exposing sub-tropical regions to prolonged periods of aridity. As leopards are highly mobile and adaptable, it is unlikely that physical barriers such a mountains and rivers impede dispersal and gene flow over the long term. On the other hand, changes in the abundance and availability of resources (such as food, water and suitable habitat) are known to decrease leopard densities and, potentially, reduce gene flow (Hinde et al., 2023). Leopard populations that inhabit arid conditions (such as the Namib desert, the Kalahari and the Arabian Peninsula (Stuart, 1975; Bothma, 1998; Jacobson et al., 2016)) tend to exist in much lower densities. Densities of leopard in arid environments range from of 0.5 leopards per 100 km² for central Namibia and 1.3–1.5 leopards per 100 km² in the Kalahari. These are significantly lower than densities observed in other parts of Africa, which can range from 2.4 leopards per 100 km² to 12.7 leopards per 100 km² (Stander et al., 1997; Balme, Slotow & Hunter, 2010; Maputla, Chimimba & Ferreira, 2013; Strampelli et al., 2018; Morris et al., 2022). Leopard home ranges, which are largely governed by vegetation productivity, prey density and the density of other leopards (Santos et al., 2019; Snider et al., 2021; le Roex et al., 2022), are also much larger in arid regions. Therefore, while increased mid-Pleistocene aridity may not have prevented gene flow among African leopards, we suggest that it must have attenuated female gene flow enough to promote the divergence of PAR-I and PAR-II lineages.

The oldest split within PAR-I was in western Africa, but in the rest of Africa, PAR-I remained relatively undifferentiated throughout the late Pleistocene. Demographic analyses suggest that PAR-I underwent a population expansion (Table 3) beginning around 200 Ka (Fig. 4). The presence of one highly frequent haplotype that gave rise to all other non-West African PAR-I haplotypes is fully in line with this hypothesis. Under a population expansion, genetic drift is less likely to sort lineages across the landscape, explaining the lack of geographic structure among non-West African PAR-I leopards.

Unlike PAR-I, PAR-II was relatively well structured geographically, with much less haplotype sharing among locations and unique sets of haplotypes occurring in Central Africa, the Cape, KwaZulu-Natal and the South African Highveld (Loskop Dam Nature Reserve). Demographic analyses indicate that the effective population size of PAR-II was likely constant over the mid-Pleistocene with a gradual decline beginning around 100 Ka (Fig. 4). A major population contraction of PAR-II within southern or South Africa is unlikely from our data, especially in light of its high genetic diversity (Table 3). Nevertheless, higher drift relative to an expanding PAR-I is more likely to have led to the evolution of location-specific PAR-II haplotypes. Higher drift in PAR-II was also reflected in increased phylogenetic structure, with central African haplotypes (DRC and Zambia) diverging first, followed by the South African Lowveld, KwaZulu-Natal and finally the Eastern and Western Cape (Fig. 1). This phylogeographic structure suggests an origin for PAR-II in Central Africa and dispersal into southern Africa during the late Pleistocene.

The differing evolutionary trajectories of PAR-I and PAR-II indicate that they must have evolved largely in isolation from each other. As explained above, the only mechanism likely to attenuate gene flow in leopards is increased aridification and associated decreases in leopard population density. The aridification of the Limpopo Basin, at the junction of the southern African countries South Africa, Botswana, Zimbabwe and Mozambique (see Fig. 3), occurred approximately from 1.0 to 0.6 Ma (Caley et al., 2018) and could provide an explanation for the isolation of the more derived South African PAR-II haplotypes (95% HPD = 0.2378–0.6857 Ma) from leopards in other parts of Africa. This period of intense aridity was also significant for hominids in Africa, with the last of the Australopithecine ape men (Paranthropus robustus) disappearing from the fossil record around this time (Owen-Smith, 2021). A highly similar pattern of mtDNA genetic structure, also thought to be driven by mid-Pleistocene aridification, is apparent among Cape and grey-footed chacma baboons (Papio ursinus), with the unique haplotypes of both subspecies found only in northern South Africa (Sithaldeen, Ackermann & Bishop, 2015). Aridification in southern Africa is also suspected to have driven mtDNA structure in giraffe (Giraffa cameloepardalis, Brown et al. (2007)) and African ground squirrels (Geosciurus inauris, Herron, Waterman & Parkinson, 2005). Distinct South African mtDNA clades were also highlighted in lion (Panthera leo, Bertola et al. (2011)) and cheetah (Acinonyx jubatus, Charruau et al. (2011)) although these were not explicitly linked to Pleistocene aridification. Nevertheless, these distinct distributions of mtDNA diversity for a range of mammal species suggests genetic discontinuity in northern South Africa, potentially driven by an arid mid-Pleistocene environment.

Finally, the occurrence of PAR-I and PAR-II, beginning in north-eastern South Africa (Kruger National Park and Lowveld locations–Manyeleti GR and Andover NR) suggests that the two lineages must have come into recent secondary contact. Contemporary gene flow between PAR-I and PAR-II was also inferred from nuclear microsatellite data across the Highveld and Lowveld regions of Mpumalanga province (Morris et al., 2023). These results suggest that this part of South Africa is one of the regions of ongoing secondary contact between Africa’s two major African leopard populations. It may also be possible that PAR-I and PAR-II come into contact in the Kalahari and Namib regions (black question marks, Fig. 3B), however, low leopard population densities in these arid zones may work to attenuate gene flow and maintain the PAR-I/PAR-II divide. Only increased sampling in these regions, together with nuclear genome-scale data can help test these ideas.

Sample provenance

While museum collections assembled during the 19^th and 20^th centuries provide an invaluable resource for evolutionary and phylogeographic studies, they are sometimes limited by incomplete or incorrect metadata. We suspect that two such specimens are contained within our data set. The first is ZMUC25, a leopard labelled “South Africa–Cape”, with no collection date. This sample has a PAR-I haplotype, which is the same as that of another sample from the same museum (ZMUC24) collected in Algeria in 1850. Given that Algeria and the Cape are on opposite ends of the continent, it is highly unlikely that both sampling localities are correct. Furthermore, during the 19^th and the first half of the 20^th century, the Cape was the only point of exit for all goods leaving South Africa, regardless of where in that vast country they originated. Therefore, given that PAR-I is not observed south of Mpumalanga province, it is unlikely that the label “the Cape” is accurate or precise enough. Similarly, the sample ZMUC3980, from the same museum again, but with the label “Burundi” and dated 1880, is almost certainly erroneous. This Burundi specimen has an identical haplotype to nine other leopards from South Africa, from the Cape and KwaZulu-Natal, in a region of the PAR-II tree that is exclusively South African. There is also the chance that the location data for both ZMUC25 and ZMUC3980 are correct, and that their non-unique NADH sequences were obtained through contamination in the laboratory.