The persistence of the critically endangered Asiatic cheetah relies upon urgent connectivity protection: a landscape genetics perspective

Decreasing genetic diversity, gene flow rates and population connectivity can increase inbreeding rates and extinction risks. Asiatic cheetah is a critically endangered mammal with large home range that suffered from extreme range reduction and population decline. Their population is now fragmented into two subpopulations. We used genetic markers to estimate genetic diversity, relatedness, minimum effective population size and gene flow, and to assess population structure. Putative corridors connecting subpopulations were inferred using connectivity models based on topography, land cover and human footprints resistance variables. Individual pairwise genetic relatedness was compared with resistance values obtained from these models and with Euclidean distances between samples to assess the most important factors shaping the current genetic structure. The estimated effective population size was extremely low (Ne = 11 to 17). Both Northern and Southern Subpopulations exhibited low genetic diversity and high relatedness. Several signatures of gene flow and movement of individuals between subpopulations were observed suggesting that inferred corridors potentially connecting subpopulations are functional. However, no traces of gene flow were observed for the latest generations, maybe due to a decrease of functional connectivity in recent years. The resistance model including all variables was best related to genetic relatedness, whereas population differentiation is mostly driven by isolation by distance. The very low estimated effective population size, decreased genetic diversity, and high relatedness of Asiatic cheetah suggests that population reinforcement, removing obstacles to connectivity and boosting prey population conservation in stepping stones are urgently needed to prevent the imminent extinction of iconic biodiversity.


Introduction
Low genetic diversity, associated with low rates of gene flow, is a major factor decreasing species' survival probability and increasing extinction risks (Frankham 2005). Habitat connectivity is critical to ensure the maintenance of gene flow, which in turn contributes to the preservation of genetic diversity and baselines population genetic structure (Vandewoestijne et al. 2008). Population structure is generally the result of several variables, such as dispersal rate, social behavior, mating system, habitat connectivity and historical factors (Epps and Keyghobadi 2015;Parreira et al. 2020). Recent anthropogenic landscape transformations have extensively changed natural habitats and their connectivity, resulting in the loss of genetic diversity in several threatened species (Ellis et al. 2021;Schlaepfer et al. 2018;Willoughby et al. 2015). This is cause of immediate concern for threatened species, and particularly for large carnivores with large home ranges, small population sizes and thus low densities (Ripple et al. 2014). In face of such events, genetic studies have been used to find the most important units for conservation, to identify and prioritize corridors that ensure population connectivity, and to propose conservation strategies (Brooks et al. 2002;Casa-Marce et al. 2013;Menotti-Raymond and O'Brien 1993).
Cheetah (Acinonyx jubatus) is a critically endangered mammal with large home range and low population density, which has suffered from extreme range reduction and population decline. As a result of two ancestral population bottlenecks, the species has very low genetic diversity in comparison with other endangered species and some felids (Dobrynin et al. 2015;Menotti-Raymond and O'Brien 1993). Asiatic cheetahs (Acinonyx jubatus venaticus) have been genetically separated from African conspecifics (Charruau et al. 2011) and have experienced more severe declines in population size and range in comparison to them (Jowkar et al. 2008;Khalatbari et al. 2018). As a result of anthropogenic habitat fragmentation and transformation, the population suffered a very strong decline over the last century, from about 400 individuals in the 1970s to less than 50 individuals in the latest assessment (Asadi 1997;Khalatbari et al. 2017). The Asiatic cheetah population is currently divided into two subpopulations located around the central deserts of Iran, Northern and Southern. A third subpopulation, Western, has gone extinct very recently (Khalatbari et al. 2017) (Fig. 1).
Recent broad-scale studies on the status and conservation of cheetahs have shown that the Asiatic population is highly inbred and displays even lower genetic diversity than other cheetah subspecies (Prost et al. 2022). However, fine-scale studies focusing on the current limited geographic distribution of the Asiatic population are needed to estimate the local genetic diversity within the subspecies. For instance, the Northern Subpopulation has the largest population size and the unique recent breeding evidence in Iran (Khalatbari et al. 2017), suggesting that regional differences in genetic diversity may occur. In addition, there is no information on whether and how populations are genetically structured, if gene flow still occurs between subpopulations, and the migration rates between them. Such information is critical for making informed management decisions and to ensure survival of the subspecies.
Recent assessments of factors affecting habitat suitability identified prey availability as the most important factor limiting the distribution of the Asiatic cheetah, followed by human-and climate-related variables (Ahmadi et al. 2017;Khalatbari et al. 2018). Indeed, over the last Zamen-e Ahoo National Park (NP); and Southern Subpopulation composed by (3) Naybandan WR, (4) Kamki Bahabad Hunting Prohibited Area (HPA), (5) Bafq Protected Area (PA), (6) Ariz HPA, (7) Dareh Anjir WR, and (8) Kalmand WR. The inset shows the study area in Iran and the green dot shows the location of the recently extinct Western Subpopulation (Khalatbari et al. 2017). Roadkill individuals (red) are from Charruau et al. (2011), wild individuals (green) indicate those that were identified from scats collected in the field and captive individuals (purple) are from confiscated cheetahs decades, increasing human population and activities have caused massive landscape transformations that resulted in declines of prey populations and their range. These changes have likely influenced the population structure and the connectivity of remaining cheetah individuals. Several studies have assessed the population connectivity of Asiatic cheetahs based on distinct methodological approaches (resistant kernel and factorial least-cost path connectivity, Moqanaki and Cushman, 2016; ensemble model and circuit theory, Ahmadi et al. 2017). These studies suggested that despite the large Euclidean distance (about 600 km) between subpopulations, the core areas inhabited by subpopulations are connected through stepping stones and corridor habitats. These putative corridors can potentially be used as dispersal routes, but additional research is needed to understand whether they indeed support the functional connectivity of cheetahs and if dispersal events effectively promote gene flow between subpopulations. Such knowledge is especially important because ca. 70% of the predicted suitable cheetah habitats are located outside protected areas (Ahmadi et al. 2017;Khalatbari et al. 2018) and also because the putative corridors linking the two remaining subpopulations are threatened by the expanding road network (Mohammadi et al. 2018). Studies evaluating vehicle collision of wildlife (Mohammadi and Kaboli 2016) have identified hotspots of cheetah potential collision in some of the putative corridors. As such, knowledge on functional connectivity is critical to mitigate possible threats to population connectivity and to promote the long-term persistence of the imperilled Asiatic population.
The main aim of this study is to assess the genetic diversity and population structure of the Asiatic cheetah and to relate the population structure with human and landscape features. In particular, we aim to answer the following questions: What is the current effective population size of the population? Which is the level of genetic diversity in each subpopulation? How are subpopulations genetically structured, and is there a genetic signature of gene flow between them? How is relatedness between individuals spatially structured? Is population connectivity mostly influenced by landscape features, or are populations mainly isolated by geographic distance? Given the larger population size of the Northern Subpopulation, we hypothesize that genetic diversity is higher there than in the Southern Subpopulation. Given the ability of cheetahs to move over large distances and the vast and open desert characteristics of the study area, only with scattered human features, we hypothesize that there is still gene flow between subpopulations and that isolation by distance is the main factor contributing to potential signatures of population structure in comparison to landscape or human-related features.

Study area
The study area covers the current range of Asiatic cheetah, including protected areas with confirmed cheetah observations in recent years, distributed in the central deserts of Iran (Khalatbari et al. 2017). Observations of cheetah are concentrated in a Northern Subpopulation composed by (1) Touran Biosphere Reserve and a (2) Miandasht Wildlife Refuge (WR) and Zamen-e Ahoo National Park (NP); and a Southern Subpopulation composed by (3) Naybandan WR, (4) Kamki Bahabad Hunting Prohibited Area (HPA), (5) Bafq Protected Area (PA), (6) Ariz HPA, (7) Dareh Anjir WR, and (8) Kalmand WR (Fig. 1). Within this area, the connectivity analyses were conducted inside a rectangular polygon covering the full extent of the distribution of samples. All following analyses were conducted in the coordinate system WGS 1984 UTM Zone 40 N.

Sample collection
Samples for the genetic study had three origins: (1) tissue samples (N = 9) collected between 2000 and 2007 from road killed specimens (Charruau et al. 2011), some of them have been collected from areas where the species is now locally extinct (e.g. Kalmand WR in Given the very low humidity in the study areas, all collected scat samples were completely dry. We stored scats in individual Ziploc bags at room temperature until being transferred to the laboratory, where they were stored at -20 °C until DNA extraction. Locations of road kill specimens were extracted from CACP reports for the period of 2004 to 2019, other cheetah observations were collected and rarefied following Khalatbari et al. (2018). The geographic location of all scat samples was recorded with a Global Positioning System device on site in the WGS84 datum.

Molecular laboratory analysis
The molecular analysis of scat samples was performed in two steps: (i) identifying the species of scat depositors, and (ii) genotyping cheetah samples for individual identification. For the first step of scat analysis, the E.Z.N.A Tissue laboratory wear, UV sterilization of all equipment before and after procedures, and cleaning laboratory surfaces with bleach between DNA extraction batches. Instruments were sterilised by ethanol flaming twice in between different samples. PCRs were carried out in another dedicated room, with similar precautions. Each PCR plate included at least one PCR negative control to monitor for potential contaminations.

Genetic parameters
Consensus genotypes over four replicates were assembled following three rules: heterozygous genotypes were accepted if the same genotype was observed in two independent PCRs; homozygous genotypes were accepted if the genotype was observed in two (tissues) or three (scats; Godinho et al. 2015) independent PCRs; only one allele was accepted (and missing data for the second allele) at a locus when observed in each of two available replicates (Appendix S1). Identical genotypes among samples were identified using Gimlet v1.3.2 (Valière 2002). The same software was used to calculate genotyping error rates, including mean allelic dropout and false allele rates across loci and samples. Finally, we estimated the probability of identical genotypes being shared by chance to assess the power of our dataset to reliably differentiate among individuals in the population (Probability of Identity, PID), and among siblings (Probability of Identity among siblings, PIDsibs) (Appendix S3). Observed heterozygosity (H o ), expected heterozygosity (H e ), and unbiased heterozygosity (uH e ) for the overall population, and for each subpopulation was calculated using GenAlEx 6.5 (Peakall and Smouse 2012). The same software was also used to test for deviations from Hardy-Weinberg equilibrium (HWE) for each locus and linkage disequilibrium between pairs of loci for each subpopulation. Allelic richness and the inbreeding coefficient (F IS ) were calculated for each subpopulation using FSTAT 2.9.3 (Goudet 2001). The same software was used to calculate genetic differentiation using the theta analogue of F ST (Weir and Cockerham 1984). Calculations of individual pairwise genetic relatedness (RL) were performed using the dyadic maximum-likelihood estimator (DyadML65) implemented in Coancestry 1.0.1.7 with 100,000 bootstraps (Wang 2011). Matrix of pairwise relatedness between individuals of each subpopulation is available in Appendix S4. The selection of relatedness estimator was based on simulation results following Pacheco et al. (2017) (see Appendix S5). Individual genotypes were further used to estimate the effective population size (N e ) of cheetah population using the singlesample bias-corrected linkage disequilibrium (LD method; Waples and Do 2010) method implemented in NeEstimator V2.1 (Do et al. 2014). Estimates were performed using 0.05, DNA Kit (Omega Bio-tek, USA) was used to extract DNA following the manufacturer's instructions. DNA extractions for scats were performed with the following modifications: manufacturer's TL buffer was replaced with Gordon's lysis buffer (0.1 M Tris-HCl, 0.1 M EDTA, 0.01 M NaCl, 1% N-lauroyl sarcosine, pH 7.5-8; Maudet et al. 2004) and half of an inhibitEX® tablet (QIAGEN, GmbH, Hilden) was added after the first lysis step and before adding the OB protease to remove potential PCR inhibitors. Scat samples were extracted in batches of 24, each batch including a negative extraction control. Following extraction, DNA was amplified using the 16S_mam primers (Deagle et al. 2009), modified to contain Illumina adaptors to allow the posterior addition of sample-specific barcodes. Details of PCR conditions and primer sequences are available in Appendix S2. Successful amplifications were then prepared for Illumina sequencing and sequences were blasted against the NCBI nt database to identify the predator species (following Library preparation and Bioinformatics pipeline sections of Khalatbari et al. 2022). For road killed samples, we received the DNA extracted from tissue following Charruau et al. (2011).
For the second step of scat analysis, all samples confirmed to be of cheetah origin were processed again for DNA extraction following the protocol of Frantz et al. (2003) after applying the Guanidinium/Thiocyanate silica method (GuSCN, Boom et al. 1990), which allows higher DNA quality and quantity. DNA extractions were performed in batches of 32 samples, each batch including a negative control. The individual ID of cheetah samples was assessed using a set of 21 microsatellites developed for Felis catus (Menotti-Raymond et al. 1999) and previously tested on cheetahs (Charruau et al. 2011;Driscoll et al. 2002;Kotze et al. 2008;Marker et al. 2008). Amplifications were performed using the QIAGEN Multiplex PCR Kit in four multiplex reactions. To minimize genotyping errors, we used a pre-amplification protocol (Smith et al. 2011) and 4 replicates per sample (Taberlet et al. 1996). PCR products were separated by size on an ABI 3130xl DNA analyzer (Applied Biosystems, Foster City, California). Microsatellite alleles were scored using the GeneScan500 LIZ size standard and GeneMapper v4.0 (Applied Biosystems, Foster City, California), and checked manually. Multiplex 2 was first used on all samples to assess the DNA quality, and samples with amplification success > 70% were selected to proceed for amplification with remaining multiplexes. For sex identification, we used ZFx gene included in Multiplex 2 (Pilgrim et al. 2005). For details on multiplex combination and thermocycling parameters see Appendix S2.
To minimise potential contaminations, DNA extractions from scat samples were carried out in a dedicated, positive air pressure laboratory (CIBIO-InBio, Vairão Campus, Portugal), following strict protocols including disposable increase the influence of known values on the interpolation values (based on their distance from the predicted points), and we used standard neighbourhood search, 10 and 15 minimum and maximum neighbours, respectively, and 8 sectors of neighbourhood search. The root-mean-square error was taken as a measure of the performance of each IDW interpolation, and models providing accurate predictions should have a root-mean-square error as small as possible (Appendix S7). The resulting individual genetic relatedness maps were then converted to raster format with 0.0089 degrees (about 1 km) of spatial resolution.
A Spatial Principal Components Analysis (SPCA) was undertaken using the continuous genetic relatedness surfaces derived in the previous step. The SPCA maps were created using the Principal Components Analysis extension of ArcGIS. For each PC map, the loading scores were used as a measure of the association between surface genetic relatedness maps and the total variance accounting for each eigenvalue was used to evaluate the level of explanation of the analysis (Appendix S8). Then, the individual PC maps were combined in a single map with the Composite Bands tool of ArcGIS and resulted in a multiple-band map displaying the individual PC maps together as a Red-Green-Blue composite.

Landscape features
Landscape features that are known to putatively influence cheetah dispersal (Cozzi et al. 2013;Pettorelli et al. 2009) were considered for building the resistance map of connectivity models. These features included: (1) Land cover variables downloaded from Tateishi et al. (2003); (2) slope derived from a digital terrain model (USGS, 2006) using the Slope function from ArcGIS; (3) human footprint features such as roads, mines, cities, villages and other human-related settlements. These features were manually digitalised from satellite imagery available at Google Earth Engine, by visual inspection of images and scanning of roads, village boundaries, city boundaries, and other human settlements, such as factories, gas stations, and airports. All polygon features were converted to raster format using Arc-GIS, with 0.0089 degrees (about 1 km) pixel size. Resistance scores were attributed to each landscape feature using expert opinions based on available knowledge on ecology and behaviour of cheetahs (Ahmadi et al. 2017;Moqanaki and Cushman 2016). For details on the attribution of resistance scores see Appendix S9.

Connectivity models
To assess connectivity between individuals we used the concept of circuit theory and CIRCUITSCAPE software 0.02 and 0.01 as the lowest allele frequency allowed. The contemporary effective population size (N e ) was also estimated using the full likelihood sibship assignment implemented in COLONY V.2.0.6.6 (Wang 2009).

Population structure
We used the Bayesian clustering method implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000) to infer the most likely number of genetic groups present among samples. We evaluated from 1 to 3 groups (K = 1 to K = 3) using the admixture model with uncorrelated allele frequencies, and without prior information on sample population membership. We performed 10 runs for each K, using a burn-in period of 50,000 and a run length of 5,000,000 iterations. Structure Harvester v.0.6.94 (Earl and vonHoldt 2012) was used to estimate the optimal number of clusters (K) through the examination of Ln p(X│K) as well as the ΔK (Evanno et al. 2005;Pritchard et al. 2000). To further identify genetic patterns, we also performed a multivariate analysis using the Discriminant Analysis of Principal Components (DAPC) (Jombart et al. 2010) implemented in the adegenet package (Jombart 2008) for R (R Core Team, 2022). This method does not require data to meet Hardy-Weinberg expectations or linkage equilibrium between loci and maximizes variation between groups while minimizing the variation within groups. To infer the number of K groups, we first used the find.clusters function ranging from 1 to 12 and the Bayesian Information Criterion (BIC). We then specified the most likely number of clusters and used the DAPC function retaining 2 principal components (Script available in Appendix S6).
Furthermore, we used individual multilocus genotypes to estimate rates of recent migration between subpopulations using the program BayesAss 3.0.4 (Wilson and Rannala 2003). The calculation was conducted with 5 × 10 7 generations and sampling every 1000 steps, and the first 5 × 10 6 generations were discarded as burn-in.

Spatial patterns in genetic relatedness
Genetic relatedness between individuals was interpolated at unsampled locations. For each individual, a continuous surface of genetic relatedness to all other individuals was created with Inverse Distance Weighting (IDW) interpolation method, implemented in ArcGIS 10.5 (ESRI, 2016) in the Geostatistical Analyst extension. The IDW function assigned genetic relatedness values to unsampled locations by weighting the average of the values available at the sampled locations, considering the inverse of the distance to each sampled location when assigning weights. The interpolation parameters were the following: power was set to 5 to identify a total of 14 individual cheetahs (3 females and 11 males) distributed by four areas, corresponding to two subpopulations. In the Northern Subpopulation, 2 females and 4 males in Touran (#1 in Fig. 1; of which one male was also identified in Naybandan) and 1 male in Miandasht (#2 in Fig. 2). In the Southern Subpopulation, 2 males in Naybandan (#3 in Fig. 1; one also identified in Touran) and 4 males and 1 female in Kamki Bahabad, Bafq, Ariz and Dareh Anjir (#4 to #7 in Fig. 1).
Amplification success from scat samples was 87.3%, and the mean rate of allele dropout and false allele across markers were 4.5% and 1.0%, respectively. The power of our dataset to accurately identify individuals is illustrated by the low values achieved for the probability of identity (PID = 1.37 × 10 − 9 ; PIDsibs = 6.68 × 10 − 5 ) achieved (Appendix S3).
Structure Harvester and STRUCTURE analysis identified two genetic groups (K = 2) as the most probable number of clusters within our dataset corresponding to Northern and Southern Subpopulations (Appendix S11 and S12). Concordantly, the DAPC analysis showed two clusters separated in PCA1, corresponding to Northern and Southern Subpopulations (Appendix S13). We then used these groups of individuals for further population genetic analyses. No deviations from Hardy-Weinberg expectations were observed and no significant LD was observed between markers (Appendix S2). Genetic diversity estimates show higher values of observed heterozygosity and allelic richness for Northern Subpopulation, and both subpopulations present higher observed than expected heterozygosity ( Table 1). The average inbreeding coefficient (F IS ) was − 0.22 and − 0.14 for Northern and Southern Subpopulations, respectively. The genetic differentiation observed between the two subpopulations was 0.115, as estimated using the theta analogue of F ST . Average pairwise relatedness between individuals was similar for both subpopulations with values of about 0.25 (range: 0.00-0.74 and 0.00-0.92 for Northern and Southern subpopulations, respectively), representing the mean relatedness expected between half-siblings or grandparents and their grandchildren. The relatedness value is lower (RL = 0.17) when individuals from both subpopulations are analysed altogether ( Table 1). All except four individuals were related to at least one other individual by over 0.5 (Appendix S4).
The LD method recovered an estimate of contemporary effective population size (N e ) of 11.7 (95% confidence interval (CI): 6.7 to 21.6) and 17.2 (95% CI: 9.9 up to 35.6) using allele frequency values of 0.05 and below 0.05, respectively. The contemporary effective population size predicted by COLONY using the full likelihood method was 11 for both random mating (95% CI = 6 to 28) and non-random mating (95% CI = 5 to 30).
(McRae and Shah 2011). When several presence points were available for one individual, the centroid of all presence points was calculated to be used as focal node. However, two very distinct localities were available in the case of M6 (distant ≈ 400 km), which forced the use of two focal nodes (one in Touran and another in Naybandan; see Appendix 1). The resistance map was built by summing resistance values of landscape features and dividing them to the maximum summed value to reach values of 1 to 100. To evaluate the effects of human activities, three connectivity models estimating landscape resistance were produced: one using all variables, one including only slope and land cover variables, and another one using only human footprint features. Then, Mantel tests ) and partial-Mantel tests (using the Euclidean distance as correcting variable) were used to compare the pairwise genetic relatedness values obtained from Coancestry with the pairwise cumulative resistance distance values estimated from each resistance map obtained from the CIRCUITSCAPE analyses to assess isolation by resistance. The same values of relatedness were also compared with Euclidean distances between individuals to assess isolation by distance. Mantel and partial-Mantel tests were performed using the vegan package implemented in the R software (R Core Team 2022). Mantel's r scores and P-values were compared to evaluate if the resistance distances were more correlated with relatedness than the Euclidean distances. This allowed evaluating the strength of correlations between pairwise relatedness with Euclidian geographic distances, and resistance models based in all variables, slope and land cover, and human footprint (Appendix S10). Finally, increasing resistance values in connectivity models were taken as an indication of decreasing putative dispersal ability of cheetah in the landscape. Then, the connectivity models based on all variables were reclassified into three categories of dispersal: unsuitable, suitable and very suitable, according to the natural breaks on data. Connectivity model was overlapped with cameratrap and run-over data to evaluate if inferred corridors spatially match the location of dispersal individuals.  N = 8). From the 138 cheetah scats, we were able to achieve a consensus genotype for 104 samples with < 7.92% missing data, which observed, such as M4 and F1 (Appendix S15). The cumulative percentage of explanation of the first three PC axes was 86.4% (Appendix S8). The SPCA map suggests spatial patterns of clustering of individual pairwise genetic relatedness that group individuals and that broadly matches the putative subpopulations (Fig. 2). Identical colours in Naybandan and Touran (marked in pink in Fig. 2) suggest the occurrence of movement between these subpopulations, as observed by the recapture of individual M6 in scat samples collected in both areas.

Genetic analysis
Estimates of contemporary gene flow infer genetic signatures of past migration in four individuals: M4 and M-PB5 exhibit signals of recent migration from Northern to Southern Subpopulations, while F3 and F-AJNarina show signals of migration from Southern to Northern Subpopulations; in all cases the migration events were estimated to have occurred two generations ago (Appendix S14).
As expected, individual relatedness maps show that geographically close individuals tend to be more related in comparison to distant ones, but a few exceptions were   (Khalatbari et al. 2017) suggests that human footprint features are mostly related to the current cheetah population structure.

Discussion
Overall, we found that the current population size of Asiatic cheetah is extremely small. The population exhibits relatively low genetic diversity, although higher in the Northern Subpopulation, and very high relatedness within subpopulations suggesting alarming inbreeding. Asiatic cheetah population is genetically structured into two subpopulations whose differentiation is driven by geographic distance. However, traces of migration and gene flow show that potential functional corridors connect the two subpopulations, although this connectivity has probably decreased in the last two generations. The genetic diversity of Asiatic cheetah is apparently lower than African populations in Botswana and Tanzania

Landscape analysis
Connectivity models derived from "all variables", "slope and land cover-only resistance variables", and "human footprint-only resistance variables" show different levels of connectivity between subpopulations (Appendix S16). The reclassified connectivity model based on all variables identified 21.1%, 37.1% and 41.8% of the study area as unsuitable, suitable and very suitable for dispersal, respectively. The distribution of pixels classified as very suitable for dispersal ensures the geographical linkage between all subpopulations (Fig. 3). Euclidean geographic distance was the most significant variable to explain the differences in relatedness (P-value of 0.001 and r score of 0.524). Among the resistance values, "all variables" resistance variables have a better combination of r score and P-value both with and without Euclidean geographic distance being used as correction factor (Appendix 10), whereas "human footprint-only" resistance variable was the following one. This population of Asiatic cheetah is estimated to be less than 50 individuals, our dataset likely provides a good representation of the population, and therefore we believe that the observed genetic structure illustrates the current situation for this population. Northern and Southern Subpopulations have been separated by a large and extensive desert, Dasht-e-Kavir (Fig. 1), probably ecologically suitable for cheetahs in historical times, but not presently due to lack of prey availability (Khalatbari et al. 2018). We hypothesize that this ecological unsuitability for sustaining a permanent population of cheetah is likely the main factor contributing to the north-south pattern.
Isolation by distance was the most important factor explaining the north-south pattern observed in comparison to resistance variables analyzed alone or together (human footprint, slope and land cover variables). "Human footprint-only" variables were included in "all variables" resistance variable, and these two resistance variables were better explaining the differences in relatedness between individuals, rather than "slope and land-cover only" resistance variables. This suggests that the growing human activities is a main obstacle hindering cheetah movements. The fact that no gene flow was observed in the last two generations together with very high relatedness values within each subpopulation suggests that the connectivity between the two populations was higher in the previous generations and may have been lost or reduced in the current time. This might be a consequence of population decline or the expansion of man-made barriers in recent years (e.g. expansion of road network) or a combination of both factors.
Our models identified several corridor habitats that allow connectivity within and among Northern and Southern Subpopulations, which has also been previously identified using different datasets (Ahmadi et al. 2017;Moqanaki and Cushman 2016;Shams-esfandabad et al. 2021). Among them, only Shams-esfandabad et al. (2021) used the current distribution of cheetah, whereas the other studies included Western Subpopulation in their analysis. Corridors predicted by Shamsesfandabad et al. (2021) and those predicted by Moqanaki and Cushman (2016) are similar to our results. Corridors predicted by Ahmadi et al. (2017), which included Kavir Subpopulation in their analysis, do not predict the Dashte Kavir as a corridor and instead are encircling it. Those corridors are probably reflecting the connectivity status of Subpopulations in the past when Kavir Subpopulations still hold cheetahs. Therefore, it is likely that the local extinction of Kavir Subpopulation has an important role in the currently reduced connectivity between Northern and Southern Subpopulations. Models of climate-change impacts for Iran suggest that lowland areas have the highest risk of exposure to extreme climates , therefore restoration of habitat in areas such as Kavir National Park should (Dalton et al. 2013;Gottelli et al. 2007), it is lower than that observed in some other endangered felids (e.g., lion and tiger; Driscoll et al. 2002;Luo et al. 2004), and it is similar to genetic diversity observed in the Iberian lynx before the population recovery (Casa-Marce et al. 2013). However, we found higher genetic diversity values for the Asiatic cheetah than those reported for the highly inbreed Asiatic lion from Gir forest (Driscoll et al. 2002) (Appendix S17). Possibly, the medium-low genetic diversity here observed for the Asiatic cheetah is a consequence of its very small population size over the last decades. The higher values of H o than H e , resulting in a negative value of F IS , could be due to sampling stochasticity and random differences in allele frequencies between very small populations. Notably, the same pattern was observed for both subpopulations, despite a possible Wahlund effect contributing to heterozygote deficiency. We speculate that this pattern may suggest that cheetahs tend to select unrelated mates (Parreira et al. 2020), likely to avoid inbreeding. This would be particularly important in a population facing a serious level of inbreeding, where only four individuals show no pairwise relatedness > 0.5 to any other individual. The very large movements observed for one male (M6; c.400 km) in this study as well as in previous studies (217 km; Farhadinia et al. 2016), may also be a mechanism of inbreeding avoidance.
The estimated contemporary effective population size for the Asiatic cheetah using two different methods varied between N e = 11 to 17, which is in line with the last assessment of population size based on direct observations and camera trapping, in which 26 individuals (adults and cubs) were reported (Khalatbari et al. 2017). The observed concordance between both estimates follows expectations, with population size and its evolutionary analogue N e suggesting an extremely low population size. As expected, the genetic diversity in the Northern Subpopulation was higher than in the Southern, likely illustrating the higher population size of the Northern Subpopulation from where several evidences of breeding have been reported (Khalatbari et al. 2017; Department of Environment and Iranian Cheetah Society unpublished data).
There is a clear and strong latitudinal pattern in the distribution of genetic variability, retrieved by both Structure and DAPC, which identified two clusters representing Northern and Southern Subpopulations, and by estimates of relatedness as a surrogate to population affinity, which also recovered a similar north-south pattern of clustering individuals. Whereas the observed pattern could be due to spatial gaps in the distribution of samples, our sampling scheme is widely representative of the current range of the species (see Khalatbari et al. 2017Khalatbari et al. , 2018. The very low sample size available here is a common shortcoming in studies addressing threatened taxa (Appendix S17). Yet, given that the identified suitable habitats and stepping stones in western margins of Dasht-e-Kavir (Ahmadi et al. 2017) and their corridors (Malakoutikhah et al. 2020), can facilitate cheetah movements and trigger the reoccupation of Kavir National Park. Still, long term population viability needs to be fully assessed, which requires presently unavailable data, such as lifetime reproductive success and mortality rates. Survival of Asiatic cheetah depends on the conservation actions taken in the next few years to protect remaining individuals and provide safe habitats and corridors for their expansion and admixture. Defining conservation measures by national experts can provide a baseline for implementation of conservation actions.
Author contributions All authors contributed to the study conception and design. Material preparation and data collection were performed by Leili Khalatbari, Hamed Abolghasemi, Ehsan Hakimi, Taher Ghadirian and Houman Jowkar. Analysis was performed by Leili Khalatbari, Raquel Godinho, Urs Breitenmoser and José Carlos Brito. The first draft of the manuscript was written by Leili Khalatbari and all authors commented on previous versions of the manuscript. All authors read and approved the final Funding Open access funding provided by FCT|FCCN (b-on). This work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 668981, by AGRIGEN-NORTE-01-0145-FEDER-000007, supported by Norte Portugal Regional Operational Programme (NORTE2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF), and by Rufford Foundation (21607-1). L.K., R.G. and J.C.B. were supported by Fundação para a Ciência e Tecnologia (PD/BD/132429/2017, 2021.00647.CEECIND and EECINST/00014/2018/CP1512/CT0001, respectively).

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright be considered in long-term management plans for cheetahs. The observation of M6 scats in Naybandan and Touran, and the previous reports of observations of the same individual in Naybandan and Yazd areas (Farhadinia et al. 2016), are an indication that cheetahs are currently using these corridors. Traces of gene flow and individual movements were observed both in recently collected samples (scat samples) and roadkill samples from about 10 years ago, indicating movement of individuals in the recent past and in the current time. The several road collisions in areas identified as likely corridors are also further evidence of their usage and that roads are impacting cheetah movements (Fig. 3a, b and  c). These corridors need effective protection, such as underpasses and fencing on both sides of the road, in areas where suitable and very suitable habitats are intersecting the roads and where collision hotspots are predicted (Mohammadi et al. 2018).

Implications for cheetah conservation
Overall, the population of Asiatic cheetah is very small, the effective population size is extremely small, genetic diversity is low and the remaining individuals are highly related to each other. Despite suitable corridors available to connect subpopulations and the suggestion of inbreeding avoidance by cheetahs, the very low effective population size and closely related individuals may hamper natural growth of population size. Therefore, and in light of these findings, we emphasize that conservation breeding programs for Asiatic cheetah using only individuals from the current population is unlikely to solve the threats that the subspecies is facing (Khalatbari et al. 2021). On the contrary, increasing population size and genetic diversity of Asiatic cheetah by translocation of African cheetahs, allowing hybridization between different subspecies, should be considered in management plans for the species (Zecherle et al. 2020), taking into consideration the intraspecific variation among different subspecies (Bertola et al. 2022). In the meantime, considering the confirmed reports of reproduction only from Touran and knowing that this area has the highest genetic diversity, we believe that conservation efforts should prioritize the protection of individuals in this area from further loss. However, given the genetic differentiation observed for the Southern Subpopulation, it is important to ensure its protection. Preserving sperm can be a potential solution for saving their genetic diversity for future population rescue. Safeguarding roads in areas that are predicted as hotspots of collision is important to decrease chances of other individuals being killed in the roads and to increase the possibility of individuals moving between subpopulations. Restoring the population of herbivores, especially the two species of gazelle (Gazella subgutturosa and G. bennettii), in the