A Chromosome-Scale Genome Assembly of the Okapi (Okapia johnstoni)

Abstract The okapi (Okapia johnstoni), or forest giraffe, is the only species in its genus and the only extant sister group of the giraffe within the family Giraffidae. The species is one of the remaining large vertebrates surrounded by mystery because of its elusive behavior as well as the armed conflicts in the region where it occurs, making it difficult to study. Deforestation puts the okapi under constant anthropogenic pressure, and it is currently listed as “Endangered” on the IUCN Red List. Here, we present the first annotated de novo okapi genome assembly based on PacBio continuous long reads, polished with short reads, and anchored into chromosome-scale scaffolds using Hi-C proximity ligation sequencing. The final assembly (TBG_Okapi_asm_v1) has a length of 2.39 Gbp, of which 98% are represented by 28 scaffolds > 3.9 Mbp. The contig N50 of 61 Mbp and scaffold N50 of 102 Mbp, together with a BUSCO score of 94.7%, and 23 412 annotated genes, underline the high quality of the assembly. This chromosome-scale genome assembly is a valuable resource for future conservation of the species and comparative genomic studies among the giraffids and other ruminants.

The okapi (Okapia johnstoni) ( Figure 1A, also known as forest giraffe, is a monotypic species in the genus Okapia and the only extant sister group to the giraffe (Giraffa spp.) within the family Giraffidae. Today, okapi only occurs in central, northern, and eastern regions of the Democratic Republic of the Congo and occasionally crosses the border to Uganda. The species is classified as "Endangered" on the IUCN Red List based on an estimated decline to one-half of its original population size during the last 24 years, which equals only about three generations (Mallon et al. 2015). This elusive species is difficult to observe in its natural habitat, dense rainforest, and most population estimates rely on dung-based surveys, as ongoing political conflicts, illegal mining, and poaching hinder more accurate population counts.
Genomic data for the okapi is still limited to being sequenced from non-invasive sampling strategies (e.g., dung samples) and single molecular markers, such as a few mitochondrial and nuclear genes. Analyses of these verify the presence of okapis south-west of the Congo River (Stanton et al. 2016) and characterize the genetic diversity across their geographic distribution (Stanton et al. 2014). There are two draft genome assemblies available for the okapi (Agaba et al. 2016;Dudchenko et al. 2017;Chen et al. 2019). However, both assemblies were constructed using short-read sequencing libraries only and are thus highly fragmented.
The development and increasing accessibility of long-read sequencing technologies and long-range scaffolding methods has made chromosome-scale genome assemblies for many animal and plant species possible (Bickhart et al. 2017;Marrano et al. 2020;Winter et al. 2020;Rhie et al. 2021), including the okapi's closest relative, the giraffe (Farré et al. 2019;Liu et al. 2021). The availability of these high-contiguity genomes enables complex comparative genomic studies and the understanding of the genetic factors underpinning unique evolutionary traits, such as hypertension, which is the norm in the exceptionally tall animal (Liu et al. 2021). Furthermore, highquality genome assemblies often serve as reference sequences to study species for which mapping of short-read sequences (e.g., from museum collections or feces) is the only source of genetic information.
Here, we report the first high-quality de novo genome assembly of the okapi based on PacBio continuous long reads, polished with short reads. Furthermore, we use existing chromatin interaction map (Hi-C) data (Dudchenko et al. 2017) to construct a chromosome-scale genome assembly as reference genome to study okapi genetic diversity, genome evolution, selection, and future comparative and population genomic analyses. This improved assembly of a genome from the sister group to the giraffe is an important resource to study the genetic basis and evolution of the unique characteristics of giraffids and will facilitate future genome-wide conservation genomic assessments of wild okapi populations.

Biological Materials
A kidney-tissue sample from a male okapi named Moyo (International Studbook Number 0385) was collected by veterinarians during post-mortem examination at the Zoo Antwerp and preserved at −80 °C until further processing. Figure 1. The okapi, Hi-C contact map, and Gene set completeness analyses. Painting of the okapi by Jon Baldur Hlidberg (www.fauna.is) (A). Proximity ligation contact map of the scaffolded assembly TBG_Okapi_asm_v1 (B). Gene set completeness analyses (BUSCO) using both the mammalia_odb10 (C) and the cetartiodactyla_odb10 gene sets (D) for the 3 okapi assemblies and the annotation of TBG_Okapi_asm_v1. The assembly with the highest BUSCO scores for both gene sets is TBG_Okapi_asm_v1.
High-molecular-weight genomic DNA (gDNA) was extracted using a standard phenol-chloroform extraction protocol (Sambrooks J, Russel DW 2001). DNA concentration and yield were measured using the Qubit dsDNA BR Assay Kit on the Qubit Fluorometer (Thermo Fisher Scientific) and DNA integrity was evaluated using the Genomic DNA ScreenTape on the Agilent 2200 TapeStation system (Agilent Technologies).

Nucleic Acid Library Preparation and Sequencing
A SMRTbell sequencing library was prepared using the SMRTbell Express Prep kit v2.0 Protocol (Pacific Biosciences-PacBio, Menlo Park, CA, USA) and sequenced on the PacBio Sequel II system on continuous long read (CLR) mode with the Sequel II Sequencing Kit 2.0 (PacBio). For preparing a Chromium sequencing library (10X Genomics, Inc., Pleasanton, CA, USA), the gDNA was size-selected using a BluePippin System (SageScience, Inc., Beverly, MA, USA) to remove DNA molecules < 30 kbp. The size-selected gDNA was then sent to SciLifeLab (Stockholm, Sweden) for library preparation and sequencing on the Illumina NovaSeq 6000 system.

Genome-Wide Heterozygosity
A folded Site Frequency Spectrum (SFS) was used to estimate the genome-wide heterozygosity of the two chromosomescale assemblies from the BAM files generated for PSMC. First, ANGSD v.0.933 (Korneliussen et al. 2014) (flag -doSaf 1) was used to estimate site allele frequencies. We set the minimum score for both mapping and base quality to 30, used the 95 th percentile of the sample's depth distribution as maximum depth cut-off, and enabled the extended Base Alignment Quality (BAQ) adjustment (flag -baq 2). The folded SFS was then generated with ANGSD's companion program realSFS (flag -fold 1) with 200 bootstrap replicates and used to calculate heterozygosity as a percentage of heterozygous sites in R v.4.1.2 (R Core Team. 2015).

Genome Sequencing and Assembly
Sequencing on the PacBio Sequel II generated 196.7 Gbp of long-read data with a mean subread length of 9.79 kbp, yielding a coverage of approximately 79-fold. For the 10X Genomics Chromium library, we received 101.6 Gbp of raw sequencing data or 677 million short reads (41-fold).
The initial supernova run to assemble the Chromium data resulted in an assembly with a total length of only 1.9 Gbp and a scaffold N50 of 47.40 kbp. Therefore, we opted for a PacBio-based de novo assembly and used the Chromium data only for polishing. The final scaffolded and gap-closed chromosome-scale assembly (TBG_Okapi_asm_v1) has a total length of 2 387 734 269 bp in 3653 scaffolds and 3728 contigs, resulting in a contig N50 of over 61 Mbp and a scaffold N50 of more than 102 Mbp (Table 1A). The 28 largest scaffolds (>3.9 Mbp), which is five or six more than the expected haploid number of chromosomes for the okapi (2n = 44-46) (Ulbrich and Schmitt 1969;Petit et al. 1994) make up 98.0% of the total assembly length ( Figure 1B). The remaining contigs are smaller than 500 kbp. The five or six scaffolds that could not be anchored into the expected 44-46 chromosome-scale scaffolds lack enough evidence in the Hi-C contact map to be placed properly either by the algorithm of the HiRise pipeline or by manual curation. One reason for this could be the use of previously published Hi-C data from a different individual, which could exhibit a different karyotype than the individual used for the new assembly.

Assembly Completeness and Quality Assessment
GenomeScope estimated the genome size of the okapi at 2 526 089 448 bp, which is 140 Mbp larger than the total assembly length of TBG_Okapi_asm_v1 but 360 Mbp smaller than ASM166083v1_HiC. The coverage-based approach used in backmap estimated the haploid genome size at 3.07 Gbp, nearly 500 Mbp larger than the k-mer based estimate. Despite the discrepancy in genome size estimates, Merqury estimated an assembly completeness for TBG_Okapi_asm_ v1 of 91.89% and a QV of 32.79 based on k-mers, corresponding to a base error rate of 0.05%.
The gene set completeness analyses in BUSCO also suggest high completeness of the assembly identifying 94.7% of complete BUSCO genes in the Mammalia and 94.3% in the Cetartiodactyla dataset, the highest BUSCO scores of the three available assemblies ( Figure 1C,D). Furthermore, long and short reads map to the assembly with a mapping rate of 95.7% and 95.5%, respectively, and no contamination is evident in the BlobPlot (Supplementary Figure S1).

Genome-Wide Heterozygosity
The genome-wide heterozygosity values for the chromosomescale assemblies from ASM166083v1_HiC and the new assembly TBG_Okapi_asm_v1 were estimated at 0.120% (sd 9.24 × 10 -7 ) and 0.173% (sd 1.22 × 10 -6 ), respectively. These new estimates are in the same range as a previous estimate (0.132%) (Brüniche-Olsen et al. 2018) derived from the same short-read data as ASM166083v1_HiC. Even though these genomes are derived from captive individuals, the current levels of heterozygosity are likely to be representative for wild okapi in that the fully resolved captive okapi pedigree suggests no recent inbreeding. For example, Moyo, the individual sampled for TBG_Okapi_asm_v1, is the offspring of a female founder individual and an unrelated 1.5-generation male.
In comparison to other species, the heterozygosity of the two okapi individuals are around three to four times higher than the heterozygosity of wild individuals of the northern giraffe Giraffa camelopardalis (0.03-0.045%), which have a census population of approximately 5000 individuals (Coimbra et al. 2021. However, these estimates were derived using a different BAQ adjustment (-baq 1) in ANGSD, which is known to underestimate the heterozygosity by approximately three to four times (Prasad et al. 2022). Estimates from Hu et al. (Hu et al. 2020) place the heterozygosity of the okapi in the same range as other endangered species such as the Sumatran Rhinoceros Dicerorhinus sumatrensis (0.130%), the Giant Panda Ailuropoda melanoleuca (0.132%), and the Western lowland gorilla Gorilla gorilla gorilla (0. 144%). Yet, estimates of heterozygosity vary greatly depending on the methods used; hence, comparisons between studies must be evaluated with caution. Furthermore, heterozygosity alone may be an inaccurate representation of a species' level of threat of extinction from genetic depletion (Teixeira and Huber 2021), thus requiring complementary assessments of mutational load, inbreeding, and population sizes.

Repeat Annotation
Repeat annotation of TBG_Okapi_asm_v1 identified a repeat content of 43.53% or 1.04 Gbp of the assembly (Table 1B).

574
Journal of Heredity, 2022, Vol. 113, No. 5 The most abundant retroelements are Long Interspersed Nuclear Elements (LINEs), which account for 31.09% of the assembly, whereas Long Terminal Repeats (LTR) elements, Short Interspersed Nuclear Elements (SINEs), and DNA transposons each only account for < 5%. The most common LINEs in TBG_Okapi_asm_v1 are RTE/Bov-B (15.27%) and L1/LINE1 (13.70%). Repeat statistics calculated for the two previously published assemblies show no major differences in repeat content compared to TBG_Okapi_asm_ v1, except of an approximately 2% lower repeat content in ASM166083v1_HiC likely attributed to the more fragmented assembly (Supplementary Table S1).

Gene Annotation
The homology-based gene prediction identified 23 412 genes in the TBG_Okapi_asm_v1 assembly with a median gene length of 7408 bp resulting in BUSCO scores of 88.8% and 90.1% for the Cetartiodactyla and Mammalia datasets, respectively ( Figure 1C,D). Of all 55 344 predicted proteins, 54 708 (98.85%) were functionally annotated by InterProScan, and in 41 684 proteins (75.31%) at least one GO term was identified. Furthermore, 53 839 proteins (97.27%) were assigned to the Swiss-Prot database.

Genome Synteny With Other Ruminants
Synteny analyses revealed high synteny between the two okapi assemblies (Figure 2A) with a few inversions and five additional scaffolds in TBG_Okapi_asm_v1. Manual curation of the Hi-C-based scaffolding yielded no evidence for the placements of the five additional scaffolds or the correct orientation of contigs responsible for the inversions in the synteny plot. As expected from cytological studies, the comparison between okapi and the Masai giraffe genome ( Figure  2B) revealed major rearrangements. In fact, not a single okapi scaffold could be found in the giraffe without being split into multiple parts or combined into the 15 much larger giraffe chromosomes. These chromosomal rearrangements, especially in the giraffe, are caused by Robertsonian-centric fusions of acrocentric chromosomes (Huang et al. 2008;Cernohorska et al. 2013;Agaba et al. 2016;Liu et al. 2021).
In contrast, the synteny between okapi and cattle or sheep was much higher ( Figure 2C,D), despite these species being more distantly related to the okapi than the giraffe. Many okapi scaffolds showed perfect synteny with scaffolds of the non-giraffid ruminants, indicating slow karyological evolution/changes.

PSMC
The estimation of effective population size (Ne) over time with PSMC revealed different trajectories for the two okapi genome assemblies TBG_okapi_asm_v1 and ASM166083v1_ HiC. For recent times (10-250 kya), both assemblies show congruent demographic patterns, with a maximum Ne around 40 kya. However, they diverge drastically further in the past (250 kya-2 Mya) reaching a 3-fold difference in Ne around 1.5 Mya (Figure 3). Differences in Ne estimates for such ancient time spans should be interpreted with caution, as they may reflect an artifact of the method. Further, both datasets are based on different library preparation methods (standard Illumina vs. 10× Genomics' Chromium), and to our knowledge, the effects of these methods on an analysis such as PSMC have not been tested. Yet, bootstrap analyses showed clear support for each of the two PSMC trajectories (Supplementary Figure S2). The third assembly (Okapi_RGD) could not be analyzed in this study, due to major differences in insert sizes in the available datasets, resulting in a much lower coverage after applying the same filtering of the BAM file as for the other datasets. However, Okapi_RGD has been previously analyzed with PSMC and showed a trajectory more similar to TBG_okapi_asm_v1 than ASM166083v1_HiC but with much more variability in the bootstrap replicates (Chen et al. 2019).

Conclusion
High-quality genome assemblies are an important resource for a variety of studies, especially for comparative genomics and conservation/ population genomics analyses. Currently, references are still missing for most non-model organisms, making the high-quality chromosome-scale genome assembly of the okapi presented here an important contribution to the field of evolutionary genomics. This is the first long-readbased assembly for this species, with the highest contiguity and completeness of BUSCO genes among the available okapi genome assemblies. Long-read-based assemblies are more likely to span complex repeat regions, such as telomeric or centromeric regions, which are usually collapsed in shortread-based assemblies, allowing for more in-depth analyses of structural variation (Gordon et al. 2016;Korlach et al. 2017;Paris et al. 2020).

Supplementary Material
Supplementary material can be found at Journal of Heredity online.

Funding
The present study was funded by the Centre for Translational Biodiversity Genomics (LOEWE-TBG) through the program "LOEWE -Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz" of Hesse's Ministry of Higher Education, Research, and the Arts.