Genome-wide diversity analyses confirm lower genetic diversity in the island species Ancient divergence between clouded leopard species predates that of leopard and lion Demographic analyses reveal distinct historical trajectories of Ne for each species

In 1821, naturalist Edward Griffith provided the first description to western science of a clouded leopard based on a skin specimen he classified as Felis nebulosa. 1 Two years later, Georges Cuvier, from a skin provided by Pierre-Me ´dard Diard from Sumatra, described a second species of clouded leopard, Felis macroscelis. 2 After examination of skulls from several different felid species in 1867, John Edward Gray assigned clouded leopards to the genus Neofelis, 3 and these were eventually subsumed into a single species, Neofelis nebulosa, a taxonomic designation that was followed for more than a century. 4,5 However, in 2006, clouded leopards were once again divided into two distinct species based on evidence from mitochondrial and nuclear DNA sequences, microsatellite and cytogenetic variation, and comparisons of pelage and skull features that all evidenced two deeply divergent and reciprocally monophyletic lineages. [6][7][8][9] These lineages were formally recognized as the mainland clouded leopard, N. nebulosa, distributed in mainland South and Southeast Asia, and the Sunda (or Diard's) clouded leopard, Neofelis diardi, found in Borneo and Sumatra. [6][7][8]10 More recently, a comprehensive evaluation of felid taxonomy affirmed this reclassification of two species and suggested that further subspecific distinctions among island populations of Sunda clouded leopards should be considered in the future. 11 Taxonomists have used biochemical and molecular genetic studies to indicate that Neofelis is sister to the Panthera lineage of cats, which includes the jaguar, lion, leopard, snow leopard, and tiger, and that these two lineages comprise the subfamily Pantherinae, the earliest diverging clade among the extant Felidae. 12,13 Previous molecular dating analyses suggest that Neofelis and Panthera diverged $5 million years ago (Mya) during the Late Miocene/Early Pliocene. [12][13][14] Because of their highly elusive natures, the ecology and habitat requirements of clouded leopards are not well known. However, clouded leopards appear to predominantly inhabit primary and secondary evergreen tropical rainforests but have also been found in other habitats such as mangrove swamps. 15,16 Like other felids, clouded leopards are hypercarnivorous, hunting a variety of terrestrial and arboreal vertebrate prey and possessing a dentition fitting such a diet, including the longest canines relative to body size among all extant felids. With long tails, broad feet and short legs that position their bodies with a low center of gravity, clouded leopards are excellent climbers, and are thought to use trees primarily as resting sites. 17 Based on analyses of communication behavior, clouded leopards appear to be solitary for most of the year except during a short mating season, similar to other solitary felid species. 17,18 Like many other terrestrial apex predators around the world, populations of the two clouded leopard species are declining and becoming fragmented throughout their respective ranges. Mainland and Sunda clouded leopards are listed as vulnerable on the IUCN Red List of Threatened Species, with mature individuals in the wild estimated to be 3,700-5,580 and 4,500, respectively. 19,20 The major threats to both species include habitat loss caused by deforestation for agricultural use, overhunting for pelts, and poaching associated with the illegal wildlife trade. 21,22 As a result of these threats, mainland clouded leopards have become locally extinct throughout their former range. As for the global ex situ population, there are 375 mainland clouded leopards housed in 108 institutions in 33 countries, especially in Europe and North America, and only $56% of the pedigree of this population is known. 23 The decline and fragmentation of clouded leopard populations raises concerns about the erosion of genetic diversity, increased probability of inbreeding, and the loss of adaptive potential. Empirical studies of genetic diversity using modern molecular genetic or genomic techniques are urgently needed for both wild and ex situ populations, but to our knowledge, no such studies have yet been completed, beyond the modeling of landscape connectivity to generate predictions of genetic diversity in Sunda clouded leopards in Sabah state, Borneo. 24 Reference genomes generated from species of conservation concern provide a first step toward assessing genome-wide genetic diversity, inbreeding, and historical demography, and facilitate the analysis of further genomic data collected from population samples. 25 Here, we present the first de novo genome assemblies for mainland and Sunda clouded leopards and evaluate their genome-wide divergence, diversity, and historical trends in effective population size.

We generated genome assemblies from one captive mainland clouded leopard (N. nebulosa) and one wild Sunda clouded leopard (N. diardi) by sequencing paired-end and mate-pair libraries to a cumulative depth of 36.7x and 34.7x, respectively (see STAR Methods). We assembled sequences de novo for both species with the MaSuRCA genome assembler. 26 For the mainland clouded leopard, we also performed Hi-C proximity ligation sequencing and assembly to generate chromosome-length scaffolds. 27 The contig and scaffold N50s for the mainland clouded leopard assembly are 76.4 kbp and 147 Mbp, respectively, and 48.9 kbp and 1.39 Mbp, respectively, for the Sunda clouded leopard assembly (see Tables S1 andS2). The contiguity of these assemblies is comparable to those of other felids including the lion, Panthera leo, (scaffold N50: 136 Mbp) 28 and the Iberian lynx, Lynx pardinus (scaffold N50: 1.52 Mbp). 29 Across 4,104 conserved mammalian genes, 30 the assemblies of the mainland and Sunda clouded leopard had Benchmarking Universal Single-Copy Orthologs (BUSCO) 31 completeness scores of 95.6 and 95.8%, respectively, indicating high gene completeness. Repetitive elements constituted 33.27% of the mainland clouded leopard genome and 32.91% of the Sunda clouded leopard genome. Homology-based gene prediction using gene annotations from human, 32 domestic dog (ROS_Cfam_1.0, NCBI Canis lupus familiaris Annotation Release 106), and domestic cat as references revealed 23,398 protein-coding genes for N. nebulosa and 23,193 for N. diardi. 33 These estimates are similar to the number of genes found in the high quality Felis_catus_9.0 assembly, which contains 24,546 genes. 33

We investigated the divergence and evolutionary history of clouded leopards within the context of felid evolution by estimating a time-scaled phylogeny using the newly generated clouded leopard genomes and those from 11 other felid species along with the spotted hyena, which was used for rooting the topologies (see STAR Methods for accession numbers and citations). We aligned 4,104 single-copy orthologs 30 (9.63 Mb) from each species and estimated phylogenies. Initial trees were estimated using maximum likelihood in IQ-TREE v.1.6.12 34 with a concatenated supermatrix (see Figure S1A) and a multispecies coalescent species tree was estimated using ASTRAL-III 35 with individual maximum likelihood gene trees as input (see Figure S1B). Both approaches yielded identical topologies within Pantherinae, with mainland and Sunda clouded leopards forming a sister clade to all Panthera species, consistent with previous studies that used smaller ($22.8-150 kb) multilocus datasets. 12,13 Phylogenetic analysis of nearly complete mitochondrial genomes, including the first such genome reported for the Sunda clouded leopard, yielded a similar deep split between Neofelis and Panthera (see Figure S2). However, relationships among species in the latter genus were discordant relative to the nuclear species tree, particularly regarding the placement of the snow leopard. Such discordance in mitochondrial phylogenies has previously been ascribed to sex-biased asymmetries in dispersal, hybrid sterility, and gene flow via historical hybridization in felids. 13 Although previous studies have noted signatures consistent with numt contamination in published mitochondrial felid genomes via segments of unusual divergence, 13 we found no such evidence of numt contamination in our newly assembled clouded leopard mitochondrial genomes using the methods outlined in Li et al. 2016. 13 Next, we used the supermatrix alignment, the species tree topology inferred with ASTRAL-III, and 11 fossilbased calibration priors and secondary priors for node ages to estimate the time to most recent common ancestor (TMRCA) among the 13 felid species in our dataset using a relaxed molecular clock implemented in MCMCTree from the PAML 4.8 package. 36 We found that mainland and Sunda clouded leopards diverged around 5.1 Mya (95% credibility interval or CI = 3.8-6.5 Mya) (Figure 2A). Although there is a wide interval of uncertainty around the point estimate, this age places the divergence of these species between the Messinian stage of the Late Miocene and the Zanclean stage of the Early Pliocene, coinciding with a highly dynamic period of sea level changes, including one of the largest marine transgressions (+50 m) at the Miocene-Pliocene boundary, $5.3 Mya. 37,38 Periods of marine transgressions would have facilitated the isolation of Sundaland from mainland Asia, assuming a model in which an ancestral population of clouded leopards was already distributed in the two regions. Our mean divergence age of 5.1 Mya is much earlier than previous estimates of 1.41 Mya and 2.86 Mya derived using smaller mitochondrial and nuclear datasets, 6,14 but is consistent with phylogenetic and historical biogeographic analyses based on combined fossil, morphological, and molecular evidence showing that a majority of pantherine lineages diversified during the Miocene and that the earliest splits among the lineages, including Neofelis, occurred in Asia. 39 The estimated divergence time between the mainland and Sunda clouded leopards is consistent with the intraspecific or interspecific ages of other forest-dependent mammal or vertebrate taxa that diversified between mainland Asia and Sundaland. 38,43 The divergence between the two clouded leopard species predates the intraspecific divergences observed in other felid species with mainland Asia and Sundaic lineages, 44 as well as the sister pair of the leopard (Panthera pardus) and lion (P. leo) at 3.5 Mya (95% CI = 2.5-4.2 Mya; Figure 1), which show more ecological and morphological differences. Despite being morphologically similar, our results suggest the split between mainland and Sunda clouded leopards is among the oldest intrageneric divergences among the extant Felidae. The split between Neofelis and Panthera is estimated at $6.8 Mya (95% CI = 5.5-8.2 Mya), which falls within the range of ages (5.5-8.3 Mya) estimated for the same node in a recent phylogenomic analysis of the extinct scimitar-toothed cat (Homotherium latidens) in relation to other extant felid species. 45

Given the deep divergence age we found between the two species of clouded leopards and because the mainland clouded leopard occupies a wider geographic range on mainland Asia relative to the insular Sunda clouded leopard, we predicted that the historical trajectory of effective population size (Ne) for the two species would be distinct and show an overall lower mean Ne for the latter species, as species might be expected to have smaller effective population sizes and lower genetic variation compared to mainland lineages. 46,47 As predicted, the Sunda clouded leopard exhibited an overall reduced Ne throughout its history relative to the mainland clouded leopard when estimated using the pairwise sequentially Markovian coalescent model (PSMC) 48 (Figure 2A). We note that the genome of the Sunda clouded leopard presented here is derived from a sample collected from a wild individual in the Malaysian state of Sabah in Borneo. Although the total estimated census size for this species is estimated to be $4,500 across both Borneo and Sumatra based on niche modeling analyses, these same analyses suggest a census size across suitable (primarily protected) habitats in Borneo of $3,800 individuals, but only $700-1,000 individuals in Sabah state, depending on the information input source used for the modeling. 24,49 Moreover, Sunda clouded leopards are thought to be patchily distributed, 24,50 which could have reduced their effective population size in more recent times.

We also tested multiple alternative mutation rates and generation lengths for comparison (see Figure S3). Our PSMC results suggest an effective population size of 3,400-18,859 for the Sunda clouded leopard and 20,048-148,718 for the mainland clouded leopard during the last 100,000 years. To compare effective population sizes for each species throughout their history, we estimated the mean effective population sizes of 13,335 and 57,255 for the Sunda and mainland clouded leopard respectively (Figure 2A). We observed a notable increase in Ne beginning around 80 Kya in the mainland clouded leopard, which is usually interpreted as a signal of population expansion. 51 However, this could also be caused by admixture among formerly isolated and structured populations. 52,53 As our mainland clouded leopard genome is derived from a zoo animal with a pedigree that includes an unknown founder history, we cannot exclude the possibility of an admixture event in this individual's history. Nonetheless, the separate Ne trajectories suggest long-term isolation between these two species.

We evaluated and visualized the density of single nucleotide variants (SNVs) across the genomes of mainland and Sunda clouded leopards to understand how genetic diversity is distributed between the two species 54 (Figure 3A). To better compare the location of SNVs in both genomes, we produced chromosome-level assemblies from the short-read assemblies of both species (see STAR Methods). We observe large blocks of low density in the Sunda clouded leopard on chromosomes A1, A2, and B3, but also several shorter regions of high density on, e.g., chromosomes A3, D3, D4, and E2. The SNV density plotted across the mainland clouded leopard was less variable (Figure 3). We estimated lower autosome-wide heterozygosity in N. diardi (0.00049) than in N. nebulosa (0.00052). Compared to other felid species for which autosome-wide heterozygosity has been estimated, clouded leopards have estimates of genetic variation comparable to that observed in cheetah (Acinonyx jubatus) and lower than most other species within the Figure 1. Phylogenetic relationships and divergence times of clouded leopards relative to other felids Time-scaled nuclear phylogeny of 13 felid species, which was rooted using the spotted hyena (Crocuta crocuta) as the outgroup (see STAR Methods). Tree topology inferred using ASTRAL-III and divergence times were estimated using MCMCTree, along with eleven secondary priors (see Table S3). Teal horizontal bars represent 95% confidence intervals. Purple bar represents the genus Neofelis. The lower figure shows changes in sea level (m) over the last 16 million years. 37 See also Figures S1 andS2.

subfamily of Pantherinae, with the exception of the snow leopard. 45 The lower heterozygosity in the Sunda clouded leopard is consistent with the lower effective population size estimated from the PSMC analyses.

The mainland clouded leopard sample was taken from a North American ex situ population, which included a founder size of 18 individuals. 55 However, gene diversity estimates of descendants based on the pedigree suggest a founder genome equivalent (FGE) as low as four individuals (FGE = 4.14). 55 Establishing compatible pairs of mainland clouded leopards of breeding age in captivity is extremely challenging, which has in some cases resulted in line breeding within the ex situ population. Consequently, the North American Species Survival Plan for mainland clouded leopards has been ''facing a demographic crisis'' because of the challenges of pairing breeding-age males and females. 56 Despite best efforts in maintaining genetic diversity in captive-bred populations, some inbreeding is still likely to occur and should be considered in this analysis. Of interest, a large chromosomal inversion was identified in the genome of another captive individual of N. nebulosa, a feature that warrants further investigation in the context of the captive breeding program (https://www.dnazoo.org/post/cloudy-with-a-chance-of-leopards).

As populations of clouded leopards continue to dwindle, understanding their individual genetic health and evolutionary and demographic histories has become crucial to their long-term sustainability. Here, we observe that mainland and Sunda clouded leopards have differing demographic histories that reveal a separation of populations and an average lower effective population size for the latter species. We also estimate a new mean divergence date of 5.1 Mya between N. nebulosa and N. diardi, which corresponds to episodes of marine transgressions that may have promoted the vicariance between populations on mainland Asia and Sundaland, as has been shown in other mammal species. The Sunda clouded leopard genome exhibits more regions of lower SNV density, which correlates with a smaller population size. In addition, this species has less average heterozygosity than observed in the mainland clouded leopard genome. These findings provide a valuable foundation for conservation managers of both captive and wild clouded leopards because they establish key differences in genetic variation between the two species. 56 Future studies should look to discover additional unique genetic differences between mainland and Sunda clouded leopards that were acquired after a 5.1 Mya divergence. It would also be beneficial to examine a larger number of individuals from both species, thereby granting a more complete picture of genetic health and demographic history.

This study presents the first chromosome-length genome assembly for the mainland clouded leopard. However, unlike the Sunda clouded leopard, which was wild caught, the mainland clouded leopard sample Cluster (SI/HPC), Smithsonian Institution (https://doi.org/10.25572/SIHPC) and the BYU Office of Research Computing high performance computing resources. We thank Hannah Landeen and Sarah Stoddard for the illustrations of the felids depicted in Figures 1 and3. We thank Christine Frandsen for helping construct the graphic abstract. The Hi-C data for Neofelis nebulosa was generaged by the DNA Zoo consortium (dnazoo.org). The DNA Zoo effort is supported by Illumina, Inc., IBM, the Pawsey Supercomputing Center, and Brigham Young University Life Sciences Inspiring Learning Internship. Unpublished genome assemblies and sequencing data for Panthera tigris, Panthera onca, Panthera uncia, Panthera pardus, Acinonyx jubatus, Crocuta crocuta, and Puma concolor are used with permission from the DNA Zoo. We also thank Jilian Fazio, the clouded leopard Species Survival Plan coordinator for the Association of Zoos and Aquariums for providing the pedigree and metadata for the mainland clouded leopard ''Ma Sing'' and sharing her insights about the ex situ population of clouded leopards in North America. We thank the Sabah Biodiversity Center for allowing us to export samples from Sunda clouded leopard individuals from Sabah (Licenses Ref JKM/MBS.1000-2/2 (391) and JKM/MNS.1000-2/3 JLD.2 (31)). We also thank Lauren Howard and the Houston Zoo for providing samples used in the Hi-C preparation of N. nebulosa. We thank three anonymous reviewers for their comments, which feel have helped us greatly improve the paper.

sample from a 3-year old wild male Sunda clouded leopard was collected into a 3.0 mL Vacutainer EDTA tube after the animal was surrendered to the Sabah Wildlife Department, Kota Kinabalu, Sabah, Malaysia underAccess License no JKM/MBS.1000-2/2 (391) and Export License no JKM/MBS.1000-2/3 JLD.2 granted by the Sabah Biodiversity Council. The sample was stored at À80 C at the Sabah Wildlife Health, Genetic and Forensic Laboratory in Kota Kinabalu until shipment to the Smithsonian National Zoological Park-Conservation Biology Institute.

We extracted genomic DNA from the two whole blood samples of the two clouded leopard species using the Qiagen DNeasy Blood & Tissue Extraction Kit (Qiagen, Valencia, California), following the manufacturer's protocols. Genomic DNA quality and quantity were evaluated with a double-stranded DNA High-Sensitivity Assay Kit (Invitrogen, Waltham, Massachusetts) and Qubit 2.0 DNA fluorometer (Life Technologies, Carlsbad, California). The genomic DNAs were then delivered to Psomagen, Inc. (Rockville, Maryland) for genomic library preparation and sequencing. Genomic DNA quantity and integrity were re-assessed using Picogreen fluorometry (Life Technologies, Carlsbad, California) and the Genomic DNA ScreenTape assay (Agilent Technologies, Santa Clara, California). Two 350-bp insert size libraries were prepared for each sample using the Illumina TruSeq DNA PCR-free kit (Illumina, San Diego, California). In addition, one 3 kb-insert size mate-pair library for each sample was prepared using the Nextera Mate Pair Library Preparation Kit (Illumina, San Diego, California). All libraries were quality checked using an Agilent Tapestation 4150 instrument and then sequenced on an Illumina HiSeq 2500 sequencer with 101 bp paired-end reads.

We trimmed the adapters from the raw read data using TrimGalore. 66 Further trimming or cleaning prior to assembling Illumina reads is not encouraged when using the MaSuRCA assembler. 26 We assembled sequencing reads from both the mainland and Sunda clouded leopard de novo with the MaSuRCA 26 assembler using the default parameters for Illumina-only assemblies. We further scaffolded the N. nebulosa assembly using Hi-C data generated by DNA Zoo to produce a chromosome-length genome assembly. A whole blood sample for in situ Hi-C preparation was donated by a female mainland clouded leopard named Rhu during a routine veterinary check-up in 2015 and provided to us by Houston Zoo. The sample was used to extract peripheral blood mononuclear cell (PBMCs) that were crosslinked and stored at À80 C. The PBCM pellet was used to prepare an in situ Hi-C library. 74 The resulting libraries (NCBI SRX7041771, SRX7041772, SRX7041774) were sequenced using the Illumina NextSeq 500 and HiSeq X Ten instruments. The Hi-C data was aligned to the draft reference assembly for N. nebulosa using Juicer 75 and assembled to chromosome-length using 3D-DNA 27 and Juicebox Assembly Tools. 76 The resulting genome assembly is available at https://www.dnazoo.org/assemblies/Neofelis_nebulosa. We assembled the mitochondrial genomes using a single paired-end library from each species using MitoFinder v1.4.1 with the-metaspades option. 77,78 We aligned the Neofelis diardi short-read assembly to the Hi-C assembly using Cactus v2019.03.01 72 and used the resulting hal alignment file as input for scaffolding with Ragout v2.3. 71 This resulted in chromosome assignments for the short-read genome and we used the Ragout-scaffolded assembly for variant calling and PSMC (see details below).

We used assembly_stats v.0.1.4 57 to generate scaffold and contig statistics for our assemblies. We then cleaned our assemblies by filtering out scaffolds below 500 bp and renaming the remaining scaffolds sequentially using Bioawk v1.0. 68 To assess the completeness of our assemblies we used Benchmark Universal Single-Copy Orthologs (BUSCO v3.0.2) to search our genomes for 4,104 mammalian orthologous genes included in the OrthoDB v.9 database. 30,31 We masked repetitive regions using Repeatmasker and the Carnivora repeat database. 58 Both genomes were annotated using GeMoMa v1.7.1, 63 which implements a homology-based gene prediction algorithm, with the human (GRCh38.p13, NCBI Homo sapiens Updated Annotation Release 109.20210514), 32 domestic dog (ROS_Cfam_1.0, NCBI Canis lupus familiaris Annotation Release 106), and domestic cat (Felis_catus_9.0, NCBI Felis catus Annotation Release 104) 33 used as references.

sampled every 10th generation until 20,000 samples of divergence time parameters were obtained, totaling 210,000 generations. We used 11 secondary priors to calibrate node ages (Table S3).

Plots showing the density of heterozygous SNVs were generated for each clouded leopard species based on the final .vcf file using VCFtools v0.1.16 67 ''snpden'' function with a window size of 1Mb using a custom script in R (https://github.com/henriquevf/snpden_plot). 73 SNV densities were then plotted onto the scaffolded assemblies of each species and then scaled to SNVs per Kbp. 54