Genomic analyses show extremely perilous conservation status of African and Asiatic cheetahs (Acinonyx jubatus)

Abstract We live in a world characterized by biodiversity loss and global environmental change. The extinction of large carnivores can have ramifying effects on ecosystems like an uncontrolled increase in wild herbivores, which in turn can have knock‐on impacts on vegetation regeneration and communities. Cheetahs (Acinonyx jubatus) serve important ecosystem functions as apex predators; yet, they are quickly heading towards an uncertain future. Threatened by habitat loss, human‐wildlife conflict and illegal trafficking, there are only approximately 7100 individuals remaining in nature. We present the most comprehensive genome‐wide analysis of cheetah phylogeography and conservation genomics to date, assembling samples from nearly the entire current and past species' range. We show that their phylogeography is more complex than previously thought, and that East African cheetahs (A. j. raineyi) are genetically distinct from Southern African individuals (A. j. jubatus), warranting their recognition as a distinct subspecies. We found strong genetic differentiation between all classically recognized subspecies, thus refuting earlier findings that cheetahs show only little differentiation. The strongest differentiation was observed between the Asiatic and all the African subspecies. We detected high inbreeding in the Critically Endangered Iranian (A. j. venaticus) and North‐western (A. j. hecki) subspecies, and show that overall cheetahs, along with snow leopards, have the lowest genome‐wide heterozygosity of all the big cats. This further emphasizes the cheetah's perilous conservation status. Our results provide novel and important information on cheetah phylogeography that can support evidence‐based conservation policy decisions to help protect this species. This is especially relevant in light of ongoing and proposed translocations across subspecies boundaries, and the increasing threats of illegal trafficking.


| INTRODUC TI ON
The current century is characterized by unprecedented global change. Climate change and human impacts disrupt the integrity of ecosystems worldwide with ramifications for faunal and floral biodiversity (Díaz et al., 2019;Eichenberg et al., 2021;Raven & Wagner, 2021). An estimated one million animal and plant species are threatened by extinction (Díaz et al., 2019). Apex predators, such as big cats, are more sensitive to environmental and climate change than species at lower trophic levels (Cheng et al., 2017;Voigt et al., 2003). Strong declines or extinctions of large carnivores can have substantial effects on ecosystems, such as an uncontrolled increase in wild herbivores, which in turn can alter plant communities and lead to shifts towards alternative ecosystem states (Beschta & Ripple, 2009) and even effect grassland net nitrogen mineralization (Frank, 2008).
The cheetah (Acinonyx jubatus, Schreber, 1775) is listed globally as a Vulnerable species by the International Union for Conservation of Nature (IUCN; Durant et al., 2015). However, as the majority of cheetahs occur outside formally protected areas, where rates of decline are likely to be elevated, Durant and colleagues argued that the cheetah meets the IUCN criteria to be categorized as Endangered (Durant et al., 2017). Furthermore, the subspecies A. j. venaticus in Iran and A. j. hecki in Northwest Africa are listed as Critically Endangered by the IUCN (Belbachir, 2008;Durant et al., 2015;Jowkar et al., 2008). Threats to cheetahs include habitat conversion and loss, persecution by pastoralists, prey declines, illegal trophy hunting, illegal trade especially as pets, and armed conflicts (Durant et al., 2014;Lindsey et al., 2011;Ray, 2005;Tricorache et al., 2018Tricorache et al., , 2021. There are approximately 7100 adult and adolescent cheetahs distributed across 33 wild subpopulations in Africa and Asia (Durant et al., 2017; Figure 1). More than half (~60%) of wild cheetahs occur in one large population in southern Africa (A. j. jubatus) (Durant et al., 2017), while A. j. venaticus is represented by fewer than 50-70 individuals and is only found in Iran (Farhadinia et al., 2017). At the end of the nineteenth century, the cheetah's distribution comprised most nonrainforest parts of Africa and much of Western and Southern Asia, from the Arabian Peninsula all the way to India, and northwards until Kazakhstan (Durant et al., 2017). However, over the past decades, the species' range declined drastically, and its current extent is probably only 9% of its historical distribution (Durant et al., 2017). Currently cheetahs are divided into four subspecies by the IUCN Cat Specialist Group (Kitchener et al., 2017), namely A. jubatus hecki (Northwest Africa), A. j. soemmeringii (Northeast Africa), A. j. jubatus (Southern and East Africa) and A. j. venaticus (Western and Southern Asia, presently found only in Iran). Krausman and Morales (2005) list A.
j. raineyi (East Africa) as a fifth subspecies, but its status is under debate (Kitchener et al., 2017). It is currently recognized as a synonym of A. j. jubatus, because of its close genetic relationship inferred from mitochondrial DNA (mtDNA; Charruau et al., 2011), a finding that has previously been reported (O'Brien et al., 1987).
Based on mtDNA and microsatellite data from up to 94 samples (OeAD) trafficking, there are only approximately 7100 individuals remaining in nature. We present the most comprehensive genome-wide analysis of cheetah phylogeography and conservation genomics to date, assembling samples from nearly the entire current and past species' range. We show that their phylogeography is more complex than previously thought, and that East African cheetahs (A. j. raineyi) are genetically distinct from Southern African individuals (A. j. jubatus), warranting their recognition as a distinct subspecies. We found strong genetic differentiation between all classically recognized subspecies, thus refuting earlier findings that cheetahs show only little differentiation. The strongest differentiation was observed between the Asiatic and all the African subspecies. We detected high inbreeding in the Critically Endangered Iranian (A. j. venaticus) and North-western (A. j. hecki) subspecies, and show that overall cheetahs, along with snow leopards, have the lowest genome-wide heterozygosity of all the big cats. This further emphasizes the cheetah's perilous conservation status.
Our results provide novel and important information on cheetah phylogeography that can support evidence-based conservation policy decisions to help protect this species. This is especially relevant in light of ongoing and proposed translocations across subspecies boundaries, and the increasing threats of illegal trafficking.

K E Y W O R D S
Acinonyx jubatus, cheetah, conservation genomics, double-digest restriction site associated DNA (ddRAD) sequencing, phylogeography including most of the cheetah's past and current range, Charruau et al. (2011) further showed that Asiatic, North-East and West African cheetahs form separate phylogenetic groups, corresponding to the currently recognized cheetah subspecies.
Recognizing the need for comprehensive conservation genomic analyses in cheetahs, we investigated genome-wide single nucleotide polymorphisms (SNPs), mtDNA and major histocompatibility complex (MHC) class II DRB immune response gene data to examine their phylogeography, and neutral and adaptive genetic diversity. We used the classical subspecies assignment based on geographic distribution as a starting point (hypothesis) to test the validity of subspecies (Sackett et al., 2014): A. j. hecki, A. j. soemmeringii, A. j. jubatus, A. j. venaticus and A. j. raineyi, acknowledging that the latter is currently not recognized by the IUCN as a separate subspecies (Kitchener et al., 2017). Discussing this information in the context of ongoing and future conservation measures, we intend that our data build the basis for comprehensive range-wide genetic monitoring of cheetahs. Moreover, it can be used to guide subspecies-specific conservation measures, as it sheds light on the genetic differentiation among cheetah F I G U R E 1 Current distribution of the five classical cheetah subspecies (after Krausman & Morales, 2005). The distribution ranges were adopted from the IUCN red list (Durant et al., 2015(Durant et al., , 2017. Subspecies were assigned to the distributions using the results of Charruau et al. (2011) and this study. Photo credits are listed in the Acknowledgements.
subspecies, and will assist with evidence-based decision-making, for example, for planned reintroduction projects and regional conservation strategies.

| MATERIAL S AND ME THODS
A detailed description of the applied methods and all Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) permits can be found in Supporting Information File 1 and   Table S1. Information about all samples used in this study can be found in Table S2.

| Data set descriptions
2.1.1 | Data set 1-nuclear DNA Data set 1 consists of 46 Acinonyx jubatus individuals used in the population genomic analyses (belonging to all five classically recognized subspecies), and 12 A. j. jubatus individuals of known parentoffspring trios only used for testing the reliability of the relatedness tests (see Table 1). Sequencing reads for three individuals in data set 1 were obtained from NCBI (Genbank: SRR2737543-SRR2737545; Dobrynin et al., 2015). The data set also includes a specimen of Puma concolor used as outgroup in the phylogenetic analyses. We extracted DNA from tissue samples of 21 modern cheetahs and of the puma using the Qiagen DNeasy Blood and Tissue kit (Qiagen) and 32 diluted blood samples using the innuPREP Blood Kit (Analytik Jena AG). We extracted the DNA of the two museum samples of A. j. hecki using the DNA extraction protocol developed by Dabney et al. (2013). All but the two museum samples were processed using the double-digest restriction-site associated DNA (ddRAD) sequencing approach (Peterson et al., 2012). The two museum specimens were sequenced for their whole genomic DNA using the method described in Meyer and Kircher (2010), as their quality was not high enough to carry out ddRAD processing. To avoid DNA contamination, we carried out all extractions in a dedicated laboratory for ancient DNA and museum samples at the University of Uppsala, Sweden.

Raw data processing
We mapped the read data against the Aci_jub_2 cheetah genome assembly (NCBI: GCA_003709585.1) using BWA-MEM (version 0.7.17-r1188) (Li & Durbin, 2010) and processed the resulting mapping files using Samtools (version 1.9) (Li et al., 2009). We subsequently assessed the mapping quality with Qualimap (Okonechnikov et al., 2016). For the two museum samples we further carried out adapter-trimming and duplicate removal using AdapterRemoval2 (Schubert et al., 2016) and Picard version 2.8.2 (https://broad insti tute.github.io/picar d/). The resulting mapping files were processed with ANGSD (Korneliussen et al., 2014), which was specifically developed for population genomic analyses of low coverage data.
We carried out filtering using SNPcleaner (version 2.24) (Fumagalli et al., 2014), by first creating a vcf file (SNP data file) for all samples using Samtools. We then filtered the SNP data for (1) the presence of no more than 25% of missing data for individual sites, (2) a maximum coverage of 120x per individual to avoid calling sites in highly repetitive regions, and (3) a minimum coverage of 3x for each individual.
The resulting sites were used as input for the downstream ANGSD analyses.

Genomic differentiation
The following analyses were carried out with the 46 individuals from data set 1 (excluding the 12 parent-offspring individuals only used for the relatedness test assessment). We carried out unsupervised principal component analyses using ANGSD and pcangsd (Meisner & Albrechtsen, 2018) and investigated population differentiation in terms of F ST with ANGSD. As some estimators of F ST can be biased by uneven sampling (see e.g., Puechmaille et al., 2016), we created three replicates for which we randomly subsampled A. j. jubatus and A. j. soemmeringii down to three individuals. All replicates resulted in highly similar estimates (Table S3), comparable to the estimates obtained from the full data set (Table S4). Furthermore, Nazareno et al. (2017) and Willing et al. (2012) showed that F ST can be reliably estimated based on small sample sizes (n > 2) when more than 1000-1500 SNPs are used in the analysis. To further investigate the reliability of our F ST estimates, we carried out a permutation test. To do so, we generated three data sets, mimicking our original sampling (four groups of three individuals and one group of two individuals), in which each group was made up of randomly selected individuals of different subspecies. We then checked the F ST distributions, the random ones and the one obtained from the original set for normality with the Shapiro-Wilk test in the R (version 3.4.3.) stats package.
As both followed a normal distribution, we next tested whether both distributions had the same variance using the var.test() function of the stats package. This showed that both distributions are independent, and we thus proceeded with an F-test (implemented in the stats package) to see whether both were significantly different from each other.
We then looked for signatures of admixture using ngsAdmix (Skotte et al., 2013). To avoid sampling biases due to different sample sizes, we subsampled A. j. jubatus and A. j. soemmeringii to three individuals to match the sample sizes of the other three subspecies.
We generated three replicates with different subsampled individuals. We further carried out a separate ngsAdmix analysis restricted to all individuals of the two subspecies A. j. soemmeringii and A. j. jubatus. We performed 50 replicates for all ngsAdmix runs ranging from k = 2 to k = 5. The results were summarized and visualized using CLUMPACK (Kopelman et al., 2015). We assessed the model fit of each k value to the data using evalAdmix (Garcia-Erill & Albrechtsen, 2020). Finally, we performed an EEMS (Petkova et al., 2016) analysis to infer effective migration rates between the African cheetah subspecies.
We carried out heterozygosity analyses for the 46 individuals using realSFS (part of the ANGSD package). Here, we did not restrict our analyses to filtered sites, to be able to compare our estimates to published genome-wide heterozygosity values for other felids (Pečnerová et al., 2021). Due to the low coverage and quality typical for degraded museum DNA, we did not include the two museum samples in this analysis.

Population structure
For the mitochondrial DNA data we aligned all sequences using Codon Code Aligner version 3.0.2 (Codon Code Corporation).
We obtained the reference sequence for the cheetah mitochondrial genome from GenBank (accession no. NC_005212.1; Burger et al., 2004). We carried out parallel analyses on the two concatenated data sets that differed in the length of the region and the number of individuals (Table S2). The first analysis comprised the larger mtDNA fragment of 929 bp amplified in 58 individuals. The second analysis was based on the shorter fragment of 681 bp and included 78 individuals from Charruau et al. (2011) and 57 examined in this study. Median-joining networks were created using Popart (Leigh & Bryant, 2015). We further investigated the molecular variance in our data using an AMOVA analysis, implemented in the software Arlequin version 3.5.2.2 (Excoffier & Lischer, 2010). Furthermore, we investigated the homology of the mtDNA regions to known NUMTs in the published cheetah genome (GCA_003709585.1) using the default settings in blastn (Altschul et al., 1990).

| Data set 3-mini-barcodes
We designed a mini-barcode approach to investigate whether all A.
j. soemmeringii carry the 3 bp deletion in the MT-ND5 gene as described in Charruau et al. (2011), and as a quick subspecies assignment test. We designed three mini-amplicons that amplify a total of 190 bp, based on diagnostic sites inferred from our mitochondrial haplotype data (data set 2). Primer sequences can be found in Table S5. We used the extracts of all A. j. soemmeringii specimens from data set 1 and sequenced the amplicons on an Illumina iSeq 100 (2 × 150 bp) following the two step sequencing protocol of Lange et al. (2014).

A. j. Soemmeringii
We mapped all sequencing reads back to the mitochondrial DNA reference for cheetahs (accession no. NC_005212.1; Burger et al., 2004) using the same approach as described for data set 1.
We then called the consensus sequences of the mapping files using ANGSD (option: -do Fasta 3), inspected them by eye, and aligned all consensus sequences using Mafft (Katoh & Standley, 2013). Then we generated median joining networks for the mini-barcode data using Popart (Leigh & Bryant, 2015).

| Data set 4-immune response genes
We sequenced the MHC class II DRB exon 2 of 46 individuals (belonging to four of the five classical subspecies; excluding A. j. hecki) to investigate their immunogenetic diversity. We used the Qiagen DNeasy Blood and Tissue kit for DNA extraction from hair and tissue samples, and the VWR PeqGold Tissue DNA Mini Kit Plus for blood samples. We carried out PCR amplifications of the target region as described in Castro-Prieto, Wachter and Sommer (2011). Indexing, multiplexing and sequencing was carried out following the Illumina Nextera XT DNA Library Prep kit reference guide. Sequencing was performed on an Illumina MiSeq (2 × 250 bp).

Adaptive immune system diversity
We mapped the reads against the MHC class II DRB exon 2 reference sequence from Castro-Prieto, Wachter and Sommer (2011) using BWA-MEM, and further processed the mapping file using Samtools.

| RE SULTS
For the ddRAD data we aimed for a minimum of 225,000 read pairs per individual, assuming the generation of 45,000 loci. We obtained 44,351 loci, and read pair counts ranged from 120,700 to 38.9 million for the different individuals, with an average read pair count of 3.58 million. For the two museum samples we aimed at 125 million read pairs and obtained 127.6 and 176 million read pairs, respectively. After filtering for coverage and missing data we retained 3743 SNPs for the analyses.

| Genomic differentiation among the five classical cheetah subspecies
Subspecies or conservation unit assignments are crucial to carry out targeted conservation efforts. Therefore, we analysed biogeographic relationships using genome-wide SNP data (3743 SNPs after filtering) for 46 individuals of the five classically recognized cheetah subspecies (see Table 1 and Table S2; data set 1). In contrast to the current and in line with the classical subspecies taxonomy, the unsupervised principal component analysis (PCA) showed five distinct genomic clusters (Figure 2a and Figure S1). These clusters correspond to the four currently recognized subspecies (A. j. jubatus, A. j. soemmeringii, A. j. hecki and A. j. venaticus) and A. j. raineyi. This is further supported by the two phylogenetic tree analyses (Figure 2b, Figure S2) and the phylogenetic network approach (Figure 2c, Figure S3). These analyses show monophyletic groups for the five subspecies (note that we could not include the two A. j. hecki samples in the ML tree analysis, which is based on genotype calls, due to their low coverage). The two phylogenetic tree analyses showed high support for these clades, with bootstraps falling in the range of 91%-100%, except for the A. j. hecki clade, which shows a bootstrap of 70% (Figure 2b and Figure S2). Overall, the two phylogenetic tree topologies were congruent with each other and placed A. j. venaticus (from Asia) as a sister group to all African subspecies. However, they differed in the placement of two individuals (ID 305 and ID306), which appeared as sisters to (i) A. j. soemmeringii using genetic distances ( Figure 2b), or (ii) A. j. raineyi in the ML tree ( Figure S2). These individuals showed mitochondrial haplotypes associated with A. j.
raineyi (see below), but clustered with A. j. soemmeringii in the PCA F I G U R E 2 Population-and phylogenomic analyses of genomewide nuclear SNP (3743) data for 46 cheetah individuals. (a) PCA analysis of all individuals of the five classical subspecies included in this study. The clustering corresponds to the morphological subspecies classification. (b) Phylogenetic relationships of representatives of the five classical cheetah subspecies using genetic distances (estimated using ngsDist and FastME; using 100 bootstraps). (c) Phylogenetic relationships of representatives for the five classical cheetah subspecies inferred by the phylogenetic network approach implemented in SplitsTree. For a fully annotated phylogenetic network see Figure S3.  NA ID432  NA ID433  BW ID430  BW ID429  ZA ID441  BW ID427  ZA ID439  BW ID428  ZA ID440  ZA ID447  ZA ID446  ZA ID445  ZA ID444  ZA ID443  ZA ID279  ZA ID438  ZA ID435  ZA ID437  ZA ID436  NA ID307  ZA ID278  BW ID442  TZ ID463  TZ ID462  TZ ID461  SD ID289  SD ID293  SO ID281  SO ID296  SO ID302  SO ID294  SO ID308  SO ID292  SO ID288  SO ID309  SO ID295  SO ID287  SO ID290  SO ID305  SO ID306  DZ ID42  EH ID75  IR ID349  IR ID348  IR  analysis based on the genome-wide SNP data. The split into five genetically distinct groups is further supported by model fit analyses based on the inferred admixture proportions (Figure 3a and Figures S4 and S5). Analyses of different subsets (at K = 5) indicated little to no signatures of admixture ( Figure 3a and Figure S4); one of the three replicate analyses suggested some, but limited, A. j.
soemmeringii ancestry in A. j. raineyi individuals and, likewise, some slight A. j. raineyi ancestry in A. j. jubatus individuals, but no signatures or patterns of admixture were consistently observed across the three different sample subsets analysed. We ran a separate admixture analysis including all individuals of A. j. soemmeringii and A.
j. jubatus, but did not detect any signatures of admixture between these two subspecies ( Figure S6). To understand patterns of gene flow among African cheetah populations, we estimated effective migration rates using EEMS (Figure 3b). This analysis indicated (1) migration between populations of the same subspecies (except for the two A. j. hecki individuals), (2)

| High inbreeding and low heterozygosity threaten the gene pool of the critically endangered Asiatic and northwest African cheetahs
We carried out inbreeding analyses using data set 1 and two dif- soemmeringii individuals (data set 1), which clearly fell within the A. j.
soemmeringii cluster based on nuclear SNP data (see Figure 2a), carried a mitochondrial haplotype belonging to the A.j. raineyi mtDNA haplogroup ( Figure S12). We found a large fraction of the total variation partitioned between the five subspecies in the hierarchical AMOVA analysis using the mtDNA data (85% between subspecies and 15% within populations).
Over all our geographical sampling we discovered four unexpected haplotypes considering their sampling locations: we found an individual from Tanzania (ID 28) and one from Zimbabwe (ID 164) with similar haplotypes to A. j. venaticus individuals (one mutation difference; Figure 4a, Figure S11). We cannot rule out that these samples were mislabelled sometime after their field collection. Also, we did not find any significant homology of these sequences to nuclear mitochondrial DNA (NUMTs) copies in the published cheetah genome (NCBI: GCA_003709585.1). Furthermore, one of the two Indian individuals included in our study (ID 425, see Table S2) shared a haplotype with an individual from Chad (ID 12, A. j. hecki, Northwest Africa; Figure 4a). This could be due to the well-documented historic cheetah trade from Africa into India (Divyabhanusinh, 1999  I b e r i a n L y n x E u r a s i a n L y n x C h e e t a h the location we would expect the Chad sample to cluster with the A. j. hecki haplogroup. More data, such as complete mitochondrial genomes or genome-wide SNP data will be necessary to determine whether these are real signals, mislabeling, or artefacts caused by the short length of the analysed mtDNA fragments.

| Adaptive immune response gene diversity
We sequenced the MHC class II DRB exon 2 of 46 individuals (belonging to four of the five classical subspecies-all but A. j. hecki; data set 4), which resulted in 13 nucleotide and nine amino acid (AA) haplotypes ( Figure 4b, Table S6). We estimated a nucleotide haplotype diversity of 0.834 (standard deviation [std]: 0.028), a nucleotide diversity (π) of 0.069 (std: 0.005), and an average of 16.2 nucleotide differences (k). The most common AA haplotype was AcjuFLA-DRB *ha16 carried by 83% of the individuals, followed by AcjuFLA-DRB *ha17 and AcjuFLA-DRB *ha20, found in 30% and 28% of the individuals, respectively (Table S6). Of the nine identified AA haplotypes, all were found in A. j. jubatus, seven in A. j. soemmeringii and five in A. j. venaticus (Table S6). AcjuFLA-DRB *ha21 and AcjuFLA-DRB *ha23 were only present in A. j. jubatus. Similar to Drake et al., 2004, we found up to four different haplotypes in a single individual (Table S6).
We further investigated whether we could expect to retrieve more haplotypes with increasing sample sizes using extrapolation (Colwell & Elsensohn, 2014; Figure S13). Even though confidence intervals were wide, especially for A. j. venaticus and A. j. soemmeringii, some trends were apparent: (1) the steady increase of haplotypes for A. j.
jubatus indicates that we have probably not saturated the sampling of MHC class II DRB exon 2 haplotypes for this subspecies; (2) the rarefaction curve of A. j. venaticus reaches a plateau at about 15 samples, which would indicate that no more than seven haplotypes are expected for this subspecies; and (3) we might be able to find more haplotypes for A. j. soemmeringii with increased number of samples ( Figure S13).

| Genomic analyses support the classical distinction of five cheetah subspecies
The understanding of subspecies structure in cheetahs has important implications for their conservation as these are often used to assign conservation units. Here, we show that our genetic results support the previously established distinction of five cheetah subspecies by Krausman and Morales (2005): A. j. jubatus, A. j. soemmeringii, A. j. venaticus, A. j. hecki and A. j. raineyi. Our reasoning follows the criteria outlined in Sackett et al. (2014). (1) Genotypic separation of putative subspecies: We found five distinct clusters in the unsupervised PCA analysis (based on nuclear SNP data), which fit the five classically recognized subspecies. Although differentiation between FST A. j. jubatus (n = 3) (n = 3) (n = 2) A. j. raineyi (n = 3)  Liu et al., 2018), lions (0.16-0.28, Smitz et al., 2018 or African leopards (0. 05-0.15, Pečnerová et al., 2021). Contrary to the mitochondrial data, A. j. raineyi displayed a closer relationship to A. j. soemmeringii than to A. j. jubatus, with which it was recently subsumed, both in the PCA and the F ST analysis using the genome-wide data (data set 1). Admixture and effective migration rate analyses indicated no or limited gene flow between the subspecies. We found two individuals assigned to A. j. soemmeringii, using genome-wide nuclear SNP data, carrying mitochondrial haplotypes of the A. j. raineyi haplogroup ( Figure S11). Interestingly, these two individuals did not show any signatures of recent admixture ( Figure S14), which could hint towards the possibility that this pattern is caused by incomplete lineage sorting rather than admixture.

A. j. hecki
More individuals are needed to formally test these different hypoth-  Figure S2). For the most part this pattern is also observed in the haplotype network reconstruction using the mitochondrial DNA data (Figure 4a; see the "Population structure using mitochondrial DNA" section of the Results for exceptions). (iv) Relatively large and significant fraction of the total variation partitioned between the subspecies: We found a large fraction of the total variation partitioned between the five subspecies in the hierarchical AMOVA analyses using the mitochondrial DNA data (data set 2; 85% between subspecies and 15% within population). We were not able to calculate hierarchical AMOVA for the genome-wide SNP data, as there are currently no available tools that can calculate this for low-coverage data.

| Phylogeography, conservation genomics and their implications
The geographic distribution of the four African cheetah subspecies is similar to that of other savannah-dwelling large African mammals.
African cheetahs are made up of Western, Northeastern, Eastern and Southern African groups. This pattern is, for example, similar to the distribution of the four currently recognized giraffe lineages (Winter et al., 2018;Coimbra et al., 2021), with a few regional differences. While the distributions of A. j. soemmeringii and A. j. raineyi meet at the border between Ethiopia, South Sudan, and Kenya, the distributions of Giraffa reticulata and G. tippelskirchi meet further south in Central Kenya. Similar to cheetahs, Southern African lions show a closer phylogenetic relationship with those in Eastern Africa than with those in Western Africa (Bertola et al., 2016).

Furthermore, population genomic analyses for lions also indicate
the presence of more than one genetic cluster in Tanzania (three for lions, and two in the case of the cheetah mitochondrial analysis). In contrast to cheetahs, where the Asian subspecies A. j. venaticus is sister to all African cheetah subspecies (Figure 2b and Figure S2), Asian lions cluster phylogenetically with North African lions (Bertola et al., 2016;de Manuel et al., 2020). Estimates of the divergence dates within cheetahs based on mitochondrial DNA or genome data vary strongly between analyses and range from 4.5 to 139 kya (see e.g., Charruau et al., 2011;Dobrynin et al., 2015;O'Brien et al., 2017;Rai et al., 2020). Also, confidence intervals in these studies often show a high degree of overlap (Rai et al., 2020) and the underlying phylogenetic trees often show limited bootstrap support for the splits separating the monophyletic subspecies clades (Charruau et al., 2011), making it difficult to decide on the sequence of subspecies divergences. Given the lack of a reliable average mutation rate estimate for cheetah nuclear DNA, the limited bootstrap support for the underlying topology even for the genome-wide SNP data, and taking the lower coverage for some of our samples into consideration, we refrained from estimating divergence times, acknowledging the limitations of our data set. High-coverage genome data from all five classically accepted subspecies will be needed to better understand the cheetah's phylogenetic history.
We found the highest levels of inbreeding in individuals of the two Critically Endangered subspecies, A. j. hecki and A. j. venaticus, further emphasizing their extremely perilous conservation status.
Moreover, A. j. venaticus also showed very low genome-wide heterozygosity. We excluded the A. j. hecki data from this analysis, because of the low DNA quality obtained from museum samples.
Of all the subspecies, A. j. jubatus showed the lowest level of inbreeding. Genome-wide heterozygosity was the highest in A. j. jubatus and A. j. raineyi. This is not surprising as A. j. jubatus makes up the largest continuous population of all cheetah subspecies and A. j. raineyi the second largest (Durant et al., 2017). We have to caution that a large fraction of our A. j. jubatus individuals (18 out of 25) originated in captivity. However, these represent purebred A. j. jubatus, and we did not observe any differences in the genome-wide heterozygosity or inbreeding estimates between our captive-bred and wild caught individuals. Our genomic analyses show that cheetahs have genome-wide heterozygosity values (0.0002-0.0005) lower than those of other endangered mid-size or large felids (data from Pečnerová et al., 2021; also noted in Dobrynin et al., 2015),

| Global change and cheetah conservation
As the cheetah serves important ecosystem functions as an apex predator, its survival is key to ensuring the healthy functioning of the ecosystems it lives in (Estes et al., 2011). The elimination or disappearance of large carnivores from fragile habitats can have knock-on effects on other faunal and floral species (Beschta & Ripple, 2009).
Threatened by global change and increasing predation pressure from humans, the protection of the last remaining cheetah populations is an imperative part of preserving natural and healthy biodiversity and habitats (Durant et al., 2017). Our findings have several implications for the conservation of cheetahs and highlight the need for further genetic studies and monitoring.

| Subspecies-specific conservation strategies
Even though our genome-wide nuclear SNP data sample set from East Africa is small (n = 3), our results indicate that A. j. jubatus and A. j. raineyi should be considered distinct at least for the purpose of management and conservation strategies (e.g., in the form of distinct conservation units). Thus, ongoing translocations of southern African cheetahs into parts of East Africa (e.g., Briers-Louw et al., 2019) should be informed by the conclusion reached in this study and may warrant closer scrutiny by conservationists. Even more so, while the phylogeographic distribution we inferred from our mtDNA data, for the most part, agrees with that described in Charruau et al. (2011), utilizing more samples, we were able to show the presence of two distinct groups within A. j. raineyi, with individuals from Tanzania and Kenya showing two different mitochondrial haplogroup assignments, falling either within the haplogroup of A. j. jubatus or forming their own haplogroup, here referred to as the A.
j. raineyi haplogroup. The three individuals from Tanzania (Dobrynin et al., 2015), which form a separate cluster using genome-wide data (data set 1), show mitochondrial haplotypes that fall within the haplogroup of A. j. jubatus. Genome-wide data for more individuals from East and Northeast Africa will be needed to better resolve the subspecies status in this region, and to better understand the complexities relating nuclear to mitochondrial DNA patterns in cheetahs.

| Development of efficient range-wide genetic monitoring strategies
We generated genome-wide data for all subspecies, which can function as a baseline for the development of reduced SNP sets, for example, for genotyping using SNP arrays or real-time PCR, enabling more cost-effective and large-scale genetic monitoring. We recently developed a SNP array for A. j. jubatus to monitor its legal and illegal trade (Magliolo et al., 2021). However, this approach was only based on individuals from one subspecies. A more generally applicable SNP typing system should include samples from all subspecies.
We caution that the presented data set should be extended by more samples of the two Critically Endangered subspecies, A. j. venaticus and A. j. hecki, and of A. j. raineyi to avoid ascertainment bias in the selected SNP sets.

| The potential for genetic monitoring of illegal trade of cheetahs
Illegal wildlife trade is a major driver of the current biodiversity loss (Maxwell et al., 2016) and affects a large portion of known plant and animal species (Fukushima et al., 2020;Scheffers et al., 2019).
Northeast Africa is a poaching hotspot for the illegal cheetah pet trade, mostly to the Gulf states (Nowell, 2014;Tricorache et al., 2018Tricorache et al., , 2021. It is also probably the region with the greatest negative impact of illegal trade on wild populations of cheetahs (Nowell, 2014).

Individuals are probably transported to the Arabian Peninsula via
Somalia and Yemen. However, the origins of these animals are poorly known. Information from interdictions and interviews with traders suggest potential origins from opportunistic collections in ethnic Somali regions such as Ethiopia and Kenya (Nowell, 2014;Tricorache et al., 2021). Interestingly, as mentioned above, Northeast Africa is the contact zone between the two subspecies A. j. soemmeringii and A. j. raineyi. Previous studies along with our current findings indicate the presence of A. j. soemmeringii in South Sudan, and the northern and central parts of Ethiopia and Somalia, and of A. j. raineyi in Kenya, Tanzania, Uganda, and the southern parts of Somalia and Ethiopia (Charruau et al., 2011;this study). Simple subspecies distinctions for illegally traded individuals and products could thus help us to quantify the respective proportion of the two subspecies in the trade, and ultimately the importance of different Northeast African countries as potential sources of origin. This can then form the basis for targeted programmes to reduce poaching and the illegal wildlife trade of cheetahs in those countries. It will also allow evidence-based decision-making for the potential release of confiscated animals into the wild, including identifying appropriate potential founders for subspecies reintroduction attempts into sites where cheetahs are currently extinct, such as Nigeria, Uganda and Rwanda. However, we caution that this aim might be complicated by possible admixture between the two subspecies or incomplete lineage sorting.

| Environmental change and immunogenetics
Immunocompetence is an important factor for the survival of a species in changing environments and is influenced by genetic factors such as MHC diversity and the environment in which individuals live (Frankham et al., 2002). For more than two decades, the cheetah has been a popular textbook example for a species with low genetic diversity, especially at MHC loci. Depleted immune gene diversity was previously supported by the cheetah's ability to accept reciprocal skin grafts from unrelated individuals (O'Brien et al., 1983) and by restriction fragment length polymorphism (RLFP) analyses of MHC class I gene transcripts and class II DRB genes (e.g., O'Brien & Yuhki, 1999). However, the findings of low MHC diversity in cheetahs has been debated after  and others (e.g., Drake et al., 2004) detected a much higher genetic diversity within MHC I loci compared to previous studies, which they attributed to a larger sample size in their study (149 cheetahs from Namibia; Castro-Prieto, Wachter & Sommer, 2011). However, they were not able to find any further MHC II-DRB haplotypes than the four previously described in Drake et al. (2004) individuals showed a higher constitutive innate immunity than sympatric leopards, which could be interpreted as a compensation for the potential lack of immunocompetence in the cheetah adaptive immune system due to their lower MHC variability (Heinrich et al., 2017). This is consistent with the lack of substantial evidence of disease events in southern and east African wild cheetah popula-

| Reintroductions of cheetahs in Asia
Several reintroduction strategies have been explored over the last years by former cheetah range countries. Frequently, the reasons for cheetah reintroductions include conservation of the species as well as expanded tourism . Ranjitsinh and Jhala (2010) identified several Indian national parks as potential candidate sites for reintroductions, although all would require extensive preparation and investment before reintroduction could be considered.
Genetic studies of regionally extinct populations would be needed to assess past genetic structure and to assign individuals to subspecies. Unfortunately, our sampling only included two individuals from India that were characterized for mtDNA only. While one individual clustered with a sample from Chad, the other one clustered with individuals assigned to A. j. venaticus (which is the suspected subspecies for cheetahs from India). Furthermore, the two Indian individuals in Charruau et al. (2011) and Rai et al. (2020) also showed A. j. venaticus haplotypes. It is well documented that imports of tamed hunting cheetahs from Northeastern (Pocock, 1939) and Eastern Africa (Divyabhanusinh, 1999) into India and the Arabian Peninsula were a regular occurrence during the European colonial era, which could explain the close relationship of one of our samples to an individual from Chad. Current proposals for cheetah restoration in India (Ranjitsinh & Jhala, 2010) suggest introducing individuals from Africa, as the last remaining A. j. venaticus representatives are highly threatened with 50-70 individuals left in the wild (Farhadinia et al., 2017). A small-scale single-locus mitochondrial genetic analysis (O'Brien et al., 2017) argued that cheetah subspecies are very closely related and that genetic distances between Asian and African cheetah subspecies are equal to those within Africa. However, our genome-wide data show that differentiation in cheetahs (average F ST of 0.38 for cheetah subspecies) is similar or even higher than that found in other large endangered felids such as tigers (0.11-0.43, Liu et al., 2018), lions (0.16-0.28, Smitz et al., 2018) and leopards (0.05-0.15, Pečnerová et al., 2021), and indicate a distinct genomewide differentiation among the African subspecies and also between Asian and African cheetahs, as we found the highest F ST values between A. j. venaticus and all the African subspecies (ranging from 0.438 to 0.497). Based on our genome-wide data and in the absence of detailed information on local and regional adaptation in different cheetah subspecies, we advise caution in releasing African cheetahs in Asia. We call for more research on the genetic and ecological differences between subspecies before resorting to introduction of African cheetah subspecies into the historical range of A. j. venaticus.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interest.

B EN EFIT-S H A R I N G S TATEM ENT
A research collaboration was developed with scientists from the countries providing genetic samples, all collaborators are included as coauthors, the results of research have been shared with the provider communities and the broader scientific community, and the research addresses a priority concern, in this case the conservation of organisms being studied. More broadly, our group is committed to international scientific partnerships, as well as institutional capacity building. All CITES permit numbers are provided in Table S1.

O PEN R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at Dryad (https://doi. org/10.5061/dryad.tx95x6b13).

DATA AVA I L A B I L I T Y S TAT E M E N T
All genomic data are deposited under the BioProject PRJNA624893 on NCBI GenBank. The mitochondrial and MHC alignments can be found on Dryad (https://doi.org/10.5061/dryad.tx95x 6b13).