Sven Winter,René Meißner,Carola Greve,Alexander Ben Hamadou,Petr Horin,Stefan Prost,Pamela A Burger
A chromosome-scale high-contiguity genome assembly of the cheetah (Acinonyx jubatus)
conservation genomics,Felidae,Hi-C,PacBio,proximity-ligation
Abstract The cheetah (Acinonyx jubatus, SCHREBER 1775) is a large felid and is considered the fastest land animal. Historically, it inhabited open grassland across Africa, the Arabian Peninsula, and southwestern Asia; however, only small and fragmented populations remain today. Here, we present a de novo genome assembly of the cheetah based on PacBio continuous long reads and Hi-C proximity ligation data. The final assembly (VMU_Ajub_asm_v1.0) has a total length of 2.38 Gb, of which 99.7% are anchored into the expected 19 chromosome-scale scaffolds. The contig and scaffold N50 values of 96.8 Mb and 144.4 Mb, respectively, a BUSCO completeness of 95.4% and a k-mer completeness of 98.4%, emphasize the high quality of the assembly. Furthermore, annotation of the assembly identified 23,622 genes and a repeat content of 40.4%. This new highly contiguous and chromosome-scale assembly will greatly benefit conservation and evolutionary genomic analyses and will be a valuable resource, e.g., to gain a detailed understanding of the function and diversity of immune response genes in felids.
A chromosome-scale high-contiguity genome assembly of the cheetah (Acinonyx jubatus)

The cheetah (Acinonyx jubatus, SCHREBER 1775) is a large felid and is considered the fastest land animal. Historically, it inhabited open grassland across Africa, the Arabian Peninsula, and southwestern Asia; however, only small and fragmented populations remain today. Here, we present a de novo genome assembly of the cheetah based on PacBio continuous long reads and Hi-C proximity ligation data. The final assembly (VMU_Ajub_asm_v1.0) has a total length of 2.38 Gb, of which 99.7% are anchored into the expected 19 chromosome-scale scaffolds. The contig and scaffold N50 values of 96.8 Mb and 144.4 Mb, respectively, a BUSCO completeness of 95.4% and a k-mer completeness of 98.4%, emphasize the high quality of the assembly. Furthermore, annotation of the assembly identified 23,622 genes and a repeat content of 40.4%. This new highly contiguous and chromosome-scale assembly will greatly benefit conservation and evolutionary genomic analyses and will be a valuable resource, e.g., to gain a detailed understanding of the function and diversity of immune response genes in felids.

The cheetah (Acinonyx jubatus, SCHREBER 1775) is a large carnivore of the cat family Felidae, in which it forms the tribe Acinonychini together with the puma (Puma concolor), and the jaguarundi (Herpailurus yagouaroundi) (Durant et al. 2021). The cheetah is known as the fastest land animal, as it reaches speeds of up to 105 km/h (Sharp 1997). Historically, it occurred in open grasslands across Africa, the Arabian Peninsula, and southwestern Asia (Durant et al. 2017). At present, it only inhabits small fractions of its former range resulting in small and fragmented populations (Durant et al. 2017). The cheetah, as a species, is currently considered “vulnerable” on the International Union for Conservation of Nature (IUCN) Red List of threatened species, with two subspecies A. j. venaticus (Iran) and A. j. hecki (Northwest Africa), being listed as “critically endangered” (Belbachir 2008; Durant et al. 2021; Farhadinia et al. 2017). The use and need for genomic analyses to support conservation decisions and management has increasingly been recognized, e.g., by the IUCN Cat Specialist Group (IUCN Cat Specialist Group 2021). Recently, conservation and evolutionary genomic analyses of the cheetah based on short-read sequences have been published (Dobrynin et al. 2015; Prost et al. 2022). However, in-depth conservation or evolutionary genomic analyses, such as the inference of runs of homozygosity (ROH) or analyses of mutational load and genetic health, greatly benefit from highly continuous reference genomes. As such, highly continuous reference genomes provide information aside from commonly used single nucleotide polymorphism or microsatellite data to evaluate a threatened species’ fitness (Wold et al. 2021) and inbreeding status (Humble et al. 2022) and play crucial roles in the development of management actions in conservation (Brandies et al. 2019). Especially ROH analyses benefit from highly continuous reference genomes. Since such a reference genome is currently unavailable for the cheetah, we sequenced and assembled a chromosome-level genome for this threatened species, with a much-improved continuity than previously reported genome assemblies.

High molecular weight genomic DNA was extracted from the blood of a 14-year-old female cheetah named “Pintada” (GAN:27869175) from Lisbon Zoo, Portugal, using the PacBio Nanobind CBB kit (PacBio, Menlo Park, CA, USA). The blood was drawn during routine veterinary procedures and immediately frozen at −20 °C. DNA concentration and molecule length were evaluated using the Qubit dsDNA BR Assay kit on the Qubit Fluorometer (Thermo Fisher Scientific) and the Genomic DNA Screen Tape on the Agilent 2200 TapeStation system (Agilent Technologies), respectively. Two sequencing libraries were prepared, one PacBio continuous long read (CLR) library using the SMRTbell Express Prep Kit 2.0 (PacBio) and one short read library using the NEBNext Ultra II FS DNA Library Kit for Illumina (New England BioLabs Inc., Ipswich, MA, USA). The long-read library was then sequenced on the PacBio Sequel IIe system in CLR mode using the Sequel II Binding Kit 2.2 (PacBio). The short-read library with an insert size of 350 bp was sequenced on the NovaSeq6000 platform (Illumina, Inc., San Diego, CA, USA), generating 150 bp paired-end reads. After receiving the data from the Sequel IIe run, we converted the PacBio subreads from BAM to FASTQ format using BAM2fastx v.1.3.0, a PacBio Secondary Analysis Tool (https://github.com/PacificBiosciences/pbbioconda; see Table 1 for a list of software used in this study). We then used Flye v.2.9 (RRID: SCR_017016) (Kolmogorov et al. 2019) to assemble the reads into a contig-level genome assembly. Flye was run with the options for raw PacBio reads and the default parameters, including one iteration of long-read polishing. The short-read data were first trimmed using fastp v.0.20.1 (RRID: SCR_016962) (Chen et al. 2018) with base correction and low complexity filter enabled to remove sequencing adaptors and polyG stretches at the end of reads. We also employed a 4 bp sliding window to detect regions of poor quality (Phred score <15). Reads were removed if they fit into one of the following categories: reads below 36 bp length, reads with >40% low-quality bases, or reads with 5 or more undetermined bases (Ns). The trimmed reads were then mapped to the assembly using bwa-mem v.0.7.17 (RRID: SCR_010910) (Li 2013). The resulting mapping file was then sorted by position in the assembly, converted to BAM format, and indexed using samtools v.1.9 (RRID: SCR_002105) (Li et al. 2009). The mapped short reads were used to further improve the base-level accuracy of the assembly with one iteration of short-read polishing with pilon v.1.23 (RRID: SCR_014731) (Walker et al. 2014). To anchor the polished contigs into chromosome-scale scaffolds, we utilized previously generated proximity ligation (Hi-C) data for the same subspecies from the DNAZoo (www.dnazoo.org, accession numbers: SRR8616936, SRR8616937). First, we followed the Arima Hi-C mapping pipeline used by the Vertebrate Genome Project (https://github.com/VGP/vgp-assembly/blob/master/pipeline/salsa/arima_mapping_pipeline.sh) for filtering and mapping of the data to the assembly. In short, the pipeline mapped the reads to the assembly using bwa-mem v.0.7.17 (RRID: SCR_010910) (Li 2013) and filtered the mapped reads with samtools v.1.14 (RRID: SCR_002105) (Li et al. 2009) based on multiple parameters such as mapping quality, read quality, and CIGAR strings. Subsequently, duplicated reads were marked and removed using the Picard v. 2.26.10 (RRID: SCR_006525) (Broad Institute 2019) tool “MarkDuplicates”. The mapped and filtered reads were then used in YaHS v.1.1 (RRID: SCR_022965) (Zhou et al. 2022) for proximity-ligation-based scaffolding. Hi-C contact maps were generated with JuicerTools v.1.22.01 (RRID: SCR_017226) (Durand et al. 2016) and used for manual curation of the scaffolded assembly in Juicebox v.1.11.08 (RRID: SCR_021172) (Durand et al. 2016). Furthermore, we run TGS-GapCloser v.1.1.1 (RRID: SCR_017633) (Xu et al. 2020) for two iterations to close gaps in the assembly and increase its contiguity. Each iteration of gap-closing utilized a random subset of approximately 25% of the long-read data to reduce computational requirements, as well as the short reads for polishing. The long-read subsets were generated from the complete dataset with seqtk v. 1.3 (RRID: SCR_018927) (Li 2018b) using the random number generator seeds 11 and 18. To evaluate the quality and completeness of the assembly, we generated assembly statistics with Quast v.5.0.2 (RRID: SCR_001228) (Gurevich et al. 2013), ran a gene set completeness analysis with BUSCO v.5.3.1 (RRID: SCR_015008) (Manni et al. 2021) using the carnivora_odb10 dataset and compared the results to the previously available chromosome-scale assembly from DNAZoo (Dobrynin et al. 2015; Dudchenko et al. 2017), which is based on an earlier draft genome (GCA_001443585.1; Dobrynin et al. 2015), and the currently best reference assembly Aci_jub_2 from GenBank (GCA_003709585.1). We also evaluated the mapping rate of both the short and long reads with QualiMap v.2.2.1 (RRID: SCR_001209) (Okonechnikov et al. 2016) after mapping the reads to our assembly with bwa-men v.0.7.17 (RRID: SCR_010910) (Li 2013) and minimap2 v.2.17 (RRID: SCR_018550) (Li 2018a), respectively. In addition, we analyzed the completeness, base-level error rate, and quality value (QV) of the assembly based on a k-mer size of 21 using Merqury v.1.1 (RRID: SCR_022964) (Rhie et al. 2020). We further evaluated potential contamination with BlobToolKit v.3.2.4 (RRID: SCR_017618) (Laetsch and Blaxter, 2017) utilizing the generated mapping files and a BLASTN v.2.11.0+ (RRID: SCR_001598) (Zhang et al. 2000) search against the NCBI Nucleotide database. In addition to the nuclear genome, we also assembled the mitochondrial genome from the short reads with GetOrganelle v1.7.5(RRID: SCR_022963) (Jin et al. 2020).

For increased accuracy during gene prediction, repeat regions in the assembly were first masked. RepeatModeler v. 2.0.1 (RRID: SCR_015027) (Flynn et al. 2020) was used to generate a de novo repeat library, which was then combined with the Felidae repeat dataset (July 2022) from RepBase (RRID: SCR_021169) (Bao et al. 2015). This custom repeat library was then used to annotate and mask the repeats in the genome using RepeatMasker v.4.1.0 (http://www.repeatmasker.org/RMDownload.html,RRID:SCR_012954). We hard-masked all interspersed repeats and soft-masked simple repeats. Genes were predicted using the homology-based gene prediction with MMseqs2 (RRID: SCR_022962) (Steinegger and Söding 2017) as an alignment tool implemented in the GeMoMa pipeline v.1.7.1 (RRID: SCR_017646) (Keilwagen et al. 2016, 2018). We used the following nine mammalian genomes and associated annotations as references: Homo sapiens (GCF_000001405.40), Mus musculus (GCF_000001635.27), Lynx canadensis (GCF_007474595.2), Canis lupus familiaris (GCF_014441545.1), Prionailuris bengalensis (GCF_016509475.1), Leopardus geoffroyi (GCF_018350155.1), Felis catus (GCF_018350175.1), Panthera tigris (GCF_018350195.1), and Panthera leo (GCF_018350215.1). Functional annotation of the predicted proteins was conducted by a BLASTP v.2.11.0+ (RRID: SCR_001010) (Zhang et al. 2000) search with an e-value cutoff of 10−6 against the Swiss-Prot database (RRID: SCR_002380; release 2021-02). Furthermore, we annotated gene ontology (GO) terms, domains, and motifs using InterProScan v.5.50.84 (RRID: SCR_005829) (Jones et al. 2014; Quevillon et al. 2005). The completeness of the predicted proteins was evaluated with BUSCO v.5.3.1 (RRID: SCR_015008) (Manni et al. 2021) using the carnivora_odb10 dataset.

We analyzed synteny between VMU_Ajub_asm_v1.0, the previously published cheetah assembly Aci_jub_2 (GCA_003709585.1), the chromosome-scale cheetah assembly from DNAZoo (Dobrynin et al. 2015; Dudchenko et al. 2017), the tiger Panthera tigris (GCF_018350195.1), two assemblies of the domestic cat Felis catus (Fca126: GCF_018350175.1, Bredemeyer et al. 2021; Fca9.1: GCA_000181335.5), the leopard cat Prionailurus bengalensis (GCA_016509475.2; Bredemeyer et al. 2021), and the Canada lynx Lynx canadensis (GCA_007474595.2; Rhie et al. 2021) using JupiterPlot v.3.8.1(RRID: SCR_022961) (Chu 2018). The closely related leopard cat was used as a reference to identify homologous cheetah chromosomes, as the chromosome structure based on G-banding is very similar to the cheetah (O'Brien et al. 2006).

Sequencing generated 341 Gb of long-read PacBio data or approximately 136-fold coverage with a mean subread length of 10,308.5 bp and approximately 18-fold (45.5 Gb) of short-read Illumina data. After assembly with Flye, polishing with pilon, proximity-ligation scaffolding with YaHS (Fig. 1A), and two iterations of gap-closing, the final assembly (VMU_Ajub_asm_v1.0) had a total length of 2.38 Gb in 198 scaffolds (including the mitochondrial genome) and a scaffold and contig N50 of 144.4 Mb and 96.8 Mb, respectively (Table 2A, Fig. 1C). The largest 19 scaffolds (>40 Mb), representing the expected haploid chromosome number of the cheetah (2n = 38) (O’Brien et al. 2006), span 99.7% of the total assembly length, resulting in a scaffold L50 of seven. VMU_Ajub_asm_v1.0 is highly contiguous and reflects a major improvement in contiguity compared with the previously available cheetah genome assembly Aci_jub_2 and the chromosome-scale one from DNA Zoo (Dobrynin et al. 2015; Dudchenko et al. 2017), as evidenced by a 569-fold and 3,007-fold larger contig NG50 (96.8 Mb vs. 170 kb vs. 32.2 kb), respectively (Table 2A). The high quality and completeness of VMU_Ajub_asm_v1.0 were also highlighted by a BUSCO completeness score of 95.4%, an increase of 0.3% and 5.7% compared with the previously available assemblies (Fig. 1B), and Merqury k-mer-based completeness of 98.4%, with an error rate of 0.013% (QV = 38.7). Furthermore, no evidence of contamination was visible in a BlobPlot (Fig. 1D).

A repeat content of 40.4% or 960.5 Mb of the sequence of VMU_Ajub_asm_v1.0 was identified by the repeat annotation (Table 2B). Long Interspersed Nuclear Elements (LINEs) were the most abundant, spanning nearly one-quarter (24.7%) of the genome, followed by Short Interspersed Nuclear Elements (SINEs) and Long Terminal Repeat (LTR) Elements with 5.45% and 5.4%, respectively. The remaining repeat classes, such as DNA Transposons, small RNA, simple repeats, etc., each accounted for less than 4%.

The homology-based gene prediction with GeMoMa identified 23,622 genes in VMU_Ajub_asm_v1.0 with a median gene length of 7,857.5 bp spanning 468.4 Mb of the total assembly length. A BUSCO score of 89.2% of identified complete Carnivora orthologous genes suggest high annotation completeness (Fig. 1B). InterProScan functionally annotated 66,775 out of the 67,405 predicted proteins (99.1%) and assigned at least one Gene Ontology (GO) term to 50,735 proteins (75.3%). In addition, more than 96.9% (65,333) of the predicted proteins were identified from the Swiss-Prot database.

JupiterPlots showed high levels of synteny between VMU_Ajub_asm_v1.0 and other felid species, as expected by the identical chromosome numbers (2n = 38) and the conserved nature of Felidae genomes (O’Brien et al. 2006) (Fig. 2). Therefore, we were able to assign chromosome names to the 19 chromosome-scale scaffolds of VMU_Ajub_asm_v1.0 according to the karyotype format commonly used for felids (O’Brien et al. 2006). We based the naming on the homologous regions of the leopard cat genome (Fig. 2C), whose chromosomes have very similar G-banding patterns as the cheetah chromosomes (O’Brien et al. 2006; Wurster-Hill and Gray 1973). Comparing VMU_Ajub_asm_v1.0 and the previous DNAZoo assembly (Fig. 2A) found structural differences only in the form of one translocation in chromosome D2 (scaffold 14) and an inversion in the smallest chromosome E4 (scaffold 19). However, both differences could potentially be scaffolding or assembly errors in one of the assemblies, despite both utilizing Hi-C data for scaffolding from the same individual. We found even fewer differences comparing VMU_Ajub_asm_v1.0 with Aci_jub_2 (Supplementary Figure 1), which was expected, as Aci_jub_2 is not chromosome-scale allowing for smaller scaffolds or contigs to be placed in syntenic positions. The most structural differences are evident between VMU_Ajub_asm_v1.0 and the most recent domestic cat assembly (Fca126, GCF_018350175.1, Fig. 2D). Yet, when compared with the previous cat assembly (Fca9.1, GCA_000181335.5, Fig. 2E), only very small rearrangements were identified, suggesting potential assembly errors in the latest cat assembly.

Highly contiguous annotated chromosome-scale genome assemblies are valuable references for evolutionary or conservation genomic analyses and enable in-depth studies on structural variation or the diversity and function of certain genes (e.g., immune response genes). However, genome assemblies from nonmodel organisms of this quality are still relatively rare. The presented new cheetah assembly VMU_Ajub_asm_v1.0, which is the first long-read-based assembly for this species, has a much-improved contiguity and will thus enable more in-depth genomic analyses for this threatened species. This genome resource provides a solid foundation to address key biological questions like understanding the process of natural selection and adaptation.

Click here for additional data file. Click here for additional data file.