Manul, also called Pallas’s cat or Rock Wildcat (Otocolobus manul), is the only representative of the genus Otocolobus, which makes up the Leopard Cat lineage along with the genus Prionailurus (Kitchener, 2017). Their habitat range is relatively wide but fragmented. Unlike other representatives of the lineage inhabiting wetland habitats, the manul lives in inhospitable cold arid montane steppes of Central Asia. Their habitat is characterized by an extreme climate with minimal rainfall and a wide range of temperatures (Ross et al., 2020). They have been recorded at altitudes up to 5,600 m (Werhahn et al., 2018). Manuls prefer rocky areas providing shelter, completely avoiding open spaces and areas with snow cover higher than 10 cm². Survival in harsh conditions is made possible by the manul’s physical characteristics. Its body structure is compact and stocky, with long shaggy fur that is double in length on its underparts - an adaptation to extremely cold winters (Ross, 2019). Short legs, low set ears and specific coat color allow the manul to blend in with its surroundings where cover options are sparse. The necessity to adapt to these conditions is probably a consequence of interspecific predation (Ross et al., 2020).

According to the IUCN Red List of endangered species, the manul is currently listed in the category “Least Concern” (Ross et al., 2020). Although their populations are probably more stable than previously thought, their current situation differs from country to country; they are classified as locally threatened in some countries and probably extinct in certain areas of its original range (Khorozyan et al., 2021). The overall population size is declining for a variety of reasons (Ross, 2019).

Individual manul populations are highly geographically fragmented, which may have negative impacts on the effective population size in terms of maintaining genomic variability. There are conservation programs in place to protect these unique cats, but information on their ecology, behavior and distribution is still rather sparse (Ross et al., 2020).

Probably due to their distinct habitat requirements, their husbandry in zoos is quite complicated and often unsuccessful despite their increasing popularity. Although diseases do not appear to be a significant threat in wild individuals, increased susceptibility to a variety of infections has been described in captive manuls (Ketz-Riley et al., 2003). Observed high mortality rates due to an acute form of toxoplasmosis appear to be a threat especially for young kittens (Brown et al., 2005; Basso et al., 2005).

The O. manul species is therefore valuable not only as a model for evolutionary and diversity studies, but it can also serve as a model for studying mechanisms of adaptation to harsh environment and mechanisms of resilience/susceptibility to diseases threatening other animal species and humans.

Whole genome sequencing (WGS) is an important tool for such studies not only in humans, but also in animals and plants (Prokop et al., 2018; Aylward et al., 2022; Pei et al., 2024). As it provides comprehensive base-by-base view of the genome, it is suitable for various applications not only in human genomics, but also in fields such as veterinary genomics, agrigenomics, ecology and biodiversity protection (Brlek et al., 2024; Wright et al., 2020). Recently, the genome of the O. manul, based on nanopore long-range sequencing of a single individual, has been published (Flack et al., 2023). To better understand inter- and intraspecific variation of the species, it is necessary to expand this information with additional resources. Using whole genome resequencing via the Illumina platform, we have obtained information on the genomes of four other manuls. It showed that all four genomes were similar in most characteristics of their variability. Information on different types of structural variants and on the extent of SNP homozygosity not available from the reference genome was retrieved. Furthermore, the four whole genome sequences, along with a recent report on unique features of the manul gene EPAS1 (Robinson et al., 2024), allowed to demonstrate how this information can be used for analyzing genes of interest, e.g., genes of adaptation.

Genomes of organisms adapted to specific environments may provide important information on their evolution and diversity, mechanisms of adaptation and/or resilience to disease. Specific adaptations to various environments were reported for various species, such as adaptation to high altitude in deer mice (Storz and Scott, 2021; Schweizer et al., 2019) or ruminants (Friedrich and Wiener, 2020), adaptation to hypoxic and cold underground environment in mole rat (Faulkes et al., 2024), longevity and unique immune adaptations in bats (Jebb et al., 2020), trypanosome tolerance (Naessens et al., 2002) or heat stress and water-deprivation stress tolerance in livestock living in tropic or arid areas (Kaushik et al., 2023; Cao et al., 2019). Information retrieved from whole genome sequences may then be used in conservation programs for endangered species as well as for optimizing environments of captive individuals (Wright et al., 2020; van der Valk et al., 2024; Brandies et al., 2019; Byrne et al., 2021). Considering the extreme climate resulting in a harsh cold and arid habitat of manuls in Central Asia, their genomes and their genetic variability is assumed to reflect the manuls’ adaptation to these specific conditions. The importance of such adaptation is illustrated by demands of keeping them in captivity (Ketz-Riley et al., 2003; Brown et al., 2005).

So far, only one whole genome sequence of O. manul has been available (OtoMan_p1.0, GCA_028564725.2). Four individual whole genome sequences obtained in this study were aligned based on this reference genome sequenced by the long-read Nanopore technology. They provide additional information on intraspecific variability, which is not available from the current reference genome assembly. As the parameters characterizing the quality of sequencing were within the standard range, the data obtained may be considered as reliable. However, the short-read sequencing technique used here has some limitations. Due to their short length, the reads may not be always successfully mapped to specific regions of the reference genome, some of them may be discarded, which may result in gaps in the sequences. The accuracy of the detection of pseudogenes, sequencing of repetitive regions, and of CNV identification is generally lower. Short-read sequences are of limited use especially in the analysis of complex genomic regions such as MHC or NKR regions, which need manual annotation (Peel et al., 2022). Despite these disadvantages, individual genomes obtained by short-read sequencing are still thoroughly explored as a useful source of additional information, especially when a long-read based reference genome is available; for example, to determine the level of diversity and signatures of selection in endangered species, which may be of use for conservation programs (Wright et al., 2020; Xia et al., 2021) or to search for loci associated with adaptation (Peng et al., 2023; Lee et al., 2018; Yu et al., 2016) or disease susceptibility (Meurs et al., 2021).

All four genomes analyzed here were similar amongst them in most characteristics. This is probably due to a sampling bias caused by limited availability of this type of samples. Tissues used for DNA extraction (manuls 4, 5 and 7) were necropsy samples originally collected for diagnostic purposes. The blood sample (manul 552) was collected for diagnostic purposes as well. The samples originated from three Czech zoos; within a zoo, the animals were related. Therefore, out of seven samples originally available, four different manul individuals (one female and three males) were selected based on available information about their origin and relatedness, with manuls 4 and 5 being littermates. Despite this primary limitation of sampling, the data obtained extend our knowledge on the overall genomic variability as well as of the structure of the individual variation represented by different types of genetic polymorphism. This information cannot be retrieved from the single reference genome OtoMan_p1.0.

We anticipated to obtain more information on the overall extent of genetic polymorphism of the species. In fact, in our four individuals sequenced, we detected a total of 3,636,571 polymorphic variants of all types on average (the mean numbers were 1,772,772 for SNPs, 1,803,852 for indels, 47,301 for SVs, and 12,647 for CNVs per individual). For comparison, 2,694,148 variants were identified in the OtoMan_p1.0 reference genome, most variants being single nucleotide variants (SNVs) and small insertions and deletions (Flack et al., 2023).

The difference in the proportion of deletions and intra-chromosomal translocations observed in manul 552 can be explained partly by its different origin, but also to a higher depth of sequencing in that individual (average depth 73.66X compared to an average of 43.38X in the others), leading to an improved detection of credibly identified variants (Sims et al., 2014). Slighter differences in the number of indels can be attributed to the varying degree of relatedness between individuals; sibling manuls 4 and 5 are more similar to each other compared to the other two.

Average whole-genome heterozygosity is another parameter, for which information from more than one reference genome allow its more accurate estimation. The average whole-genome heterozygosities of both SNP/indels and SNPs alone were almost identical between our samples and the O. manul reference genome. The average number of heterozygous SNPs per 1 Mb was slightly higher in the reference genome compared to the average of the newly sequenced manuls. Altogether, this limited dataset indicated moderate diversity of the genomes analyzed. Our estimate of whole-genome heterozygosity is limited not only by small sample sizes, but potentially also by other biological and non-biological factors, such as population structure or settings of bioinformatic analyses (rare allele filtering), respectively, that must be taken into account (Schmidt et al., 2021). Nevertheless, it does not contradict data on genetic diversity determined by genotyping 29 microsatellite loci on six wild manuls and ten manuls from a zoo-managed population, the former group showing slightly more genetic diversity (Robinson et al., 2024). The average SNP heterozygosity was 0.64 in our newly sequenced samples and the average heterozygosity of 29 microsatellite markers in wild and captive manuls was 0.54 and 0.59, respectively. Our small and partly related zoo-managed group is supposed to represent only a limited extent of existing variability both in captive and wild animals. On the other hand, taken together with the microsatellite data suffering from similar limitations (Robinson et al., 2024), the so far available data suggest that estimated nuclear heterozygosity levels of manuls as species are somewhere in the middle of the scale of mammalian species compared by Westbury et al. (Westbury et al., 2018), close to the panda.

However, at this stage of our knowledge, no conclusion on the diversity and genome-based conservation status of the manul captive and wild populations can be made, which emphasizes the need for extending genomic analyses of this understudied species. In a broader context, it should be kept in mind that the degree of vulnerability does not always correspond to the average heterozygosity of the genome (Westbury et al., 2018). Information on genome variability is thus important, but not sufficient on its own to decide on the IUCN red-list status and further management of a species (Schmidt et al., 2023). With the in-depth analysis of the four new manul genomes made available here, we add to the genomic resources of the species, which will be useful for future metapopulation studies about genome-wide diversity, genetic load and (captive) population management. Based on the IUCN Cat Specialist Group special issue on Pallas’s cat conservation (Barclay et al., 2019), manuls have not yet been integrated into any formal national action plan within their range countries, although strategic plans on the sustainable management exist. As such, more genomic information directly adds to the greater international visibility and credibility of conservation efforts in manuls as demanded by the authors. This will increase the potential to implement manul management actions into the National Biological Strategies and Action Plans that have been developed by 14 of the 16 range countries as an obligation to the Conservation on Biology Diversity (CBD 2019).

The nonsynonymous serine to proline substitution in exon 15, residue 806, within the hypoxia-related EPAS1 (endothelial PAS domain protein 1) gene was found in manuls, but not in the genomes of domestic cats or snow leopards. This is a first example of a candidate gene that could be studied in the context of manuls’ adaptation, in this particular case to hypoxic environments (Robinson et al., 2024). Our data confirmed the presence of the CCC→TCC substitution in all manuls analyzed, and its absence from the genomes of the jungle cat and black-footed cat. A BLAST search of other felid NCBI genomes did not reveal it in any other felid genome suggesting that it is really unique for this species. The study by Robinson et al. (2024) focused on variation of two exons. Our analysis covered the complete CDS and identified three additional non-synonymous substitutions unique for manuls: N84H in exon 3, S446G in exon 10 and G611S in exon 13. When we compared whole EPAS1 genomic sequences retrieved from genomes of the two felid species living in extreme but different habitats (Gray, 2014; Sliwa et al., 2013), the first part (24 bp) of the sequence annotated in the domestic cat reference genome was missing in manuls as well as in F. chaus and F. nigripes. This was then confirmed, and further variation of this part of the EPAS1 genomic sequence was observed in some of the available annotated felid sequences at NCBI. Most available felid species carry a complete EPAS1 sequence similar to the one annotated in the F. catus reference genome. Besides the three species analyzed here, variation in the length of first part of the sequence was also observed in Panthera uncia and Puma concolor. Considering the quality of annotation of most of the genomes analyzed, it is probably not an error in the annotations. Perhaps alternative splicing could be an explanation of these findings (Verta and Jacobs, 2022).

We found no polymorphic sites within the manul EPAS1 coding sequences, suggesting possible effects of purifying selection (Nielsen, 2005). Although according to SIFT, the species-specific CCC→TCC substitution should not affect the EPAS1 protein function, another non-synonymous substitution leading to interspecific differences and potentially affecting function (N84H) was identified in the manul EPAS1 gene, but not found in the two other species analyzed. Nevertheless, the importance of the interspecific variation and of the 5′UTR and 3′UTR within species variation of the EPAS1 gene remain unknown. Mutations in EPAS1 are also known in other species adapted to hypoxic environments (Lee, 2024). In human populations living in high altitudes, a reduced expression of hypoxia-inducible factor (HIF) genes protecting the body from negative effects of prolonged hypoxia, more precisely from negative effects of the mechanisms activated during hypoxia, was observed (Bigham and Lee, 2014). Interestingly, different mutations resulting in the same phenotype were observed among human populations. A mutation in the protein coding sequence was observed in the Andean population, whereas a non-coding regulatory element SNPs were identified in the Tibetan population (Lawrence et al., 2024). These findings suggest that in manuls, the non-synonymous substitutions could play a similar role in the adaptation to the high-altitude environment. In this study, we observed inter-specific and within-species variation in 5′, 3′ as well as in intronic sequences of EPAS1. The unique variants found in manuls’ lncRNA could have a significant regulatory function, similar to the situation found in the Tibetan population. At this stage however, we may only speculate about their functional importance. Expression studies are necessary to address these issues.

Altogether, this study brings the first insight into the within-species genomic variability of manuls and shows that individual whole genome sequences may serve as novel resources for studying mechanisms of adaptation and other special phenotypes characteristic for the species, which was demonstrated on an example of EPAS1 gene and lncRNA. Although the primary aim of this study was to identify genomic sequences contributing to within-species variation, we reflected first data on individual variation of both wild and captive manuls, including a newly identified non-synonymous substitution in the EPAS1 gene (Chen et al., 2009) to illustrate the potential of exploring genomes of specifically adapted species and of their comparisons with other species. The study thus underlines the importance of generating further genomic data on this understudied but in many aspects interesting species.