Reference genome of the leopard seal (Hydrurga leptonyx), a Southern Ocean apex predator

Apex predators play a key role in maintaining ecosystem structure and functions (Ripple et al., 2014, Estes et al., 2016, Enquist et al., 2020). They exert top-down control on food webs, affecting nutrient and carbon cycles, modifying habitats, and regulating the spread of disease and invasive species (reviewed in Hammerschlag et al., 2019). Hence, knowledge of the basic biology of apex predator species facilitates the prediction of their ecosystem impacts. However, it can be difficult to gather this knowledge for species that occupy remote habitats and have solitary lifestyles (Olivier et al., 2022).The leopard seal (Hydrurga leptonyx) is an apex predator primarily found at low density (Southwell et al., 2008) in the Southern Ocean and subantarctic regions (Rogers, 2009;van der Linde et al., 2022;Borras-Chavez et al., 2024). Its diverse diet includes krill, cephalopods, fish and warm-blooded prey (e.g. seabirds and other seals; Siniff and Stone, 1985;Rogers, 2009), impacting different prey populations. The resulting predation pressure can be disproportionately high, which can significantly contribute to prey population collapse (Boveng et al., 1998;Schwarz et al., 2013;Krause et al., 2022). In the last decades leopard seals have been studied more intensively, especially in subantarctic areas, and more information about their foraging and movement patterns, physiology, morphometrics, and breeding behavior has been gathered (e.g., Rogers, 2017;Staniland et al., 2018;Krause et al., 2020;Kienle et al., 2022;Sperou et al., 2023;Kienle et al., 2024). Nevertheless, the relatively scarce information available for this species constrains our understanding of its role within the Antarctic ecosystem.Genomic approaches can provide new insights into marine mammal ecology and evolution, such as understanding their spatiotemporal occurrence and abundance (Boyse et al., 2024), exploring their predator-prey dynamics (Visser et al., 2021), and investigating their demographic histories and population dynamics (Peart et al., 2020;Nebenführ et al., 2024). These methods are also growing in popularity due to the high-resolution they provide, leading to better-resolved results (Waldvogel et al., 2020). The affordability and accuracy of whole genome sequencing have advanced significantly with the development of technologies such as Illumina 10X, Oxford Nanopore, and PacBio HiFi, which are complemented by the ongoing improvement of bioinformatic tools. The use of whole genome sequencing also avoids common biases that are often encountered with specific genetic markers or genes. Studies aiming to investigate genetic variability within species or populations, or to characterize interspecific phylogenetic relationships, typically rely on reference genomes to facilitate the mapping of sequencing reads from multiple individuals (Fuentes-Pardo and Ruzzante, 2017; Theissinger et al., 2023). Reference genomes are also required for comparative genomics, transcriptomics, and epigenomics -scientific approaches that investigate species adaptation (Khudyakov et al., 2015;Yuan et al., 2021) and responses to changing environments (reviewed in Bernatchez, et al., 2024).Nuclear reference genomes have been generated for only 17 pinniped species, representing approximately half of the taxonomic diversity within this group (Supplementary Table S1). Within the Monachinae subfamily, genome assemblies are currently available for four species: the Weddell seal (Leptonychotes weddellii), Northern elephant seal (Mirounga angustirostris), Southern elephant seal (Mirounga leonina), and Hawaiian monk seal (Neomonachus schauinslandi; Supplementary Table S1). The Weddell seal is the only species from the Lobodontini tribe with a sequenced nuclear genome. This phylogenetic clade also includes the crabeater seal (Lobodon carcinophaga), Ross seal (Ommatophoca rossii), and leopard seal, all of which lack nuclear genome assemblies. Genomic resources for the leopard seal are currently limited to a mitochondrial genome (Arnason et al., 2006). Most genetic studies involving this species have employed classical genetic markers and targeted a broader context. Mitochondrial genes such as cytochrome b (cyt b), NADH dehydrogenase subunit 3 (ND3), and nuclear markers like the recombinase activating protein 1 (RAG1) have been sequenced in association with phylogenetic studies of Pinnipedia (Davis et al., 2004;Arnason et al., 2006;Fulton and Strobeck 2010). Only three studies focusing on leopard seal genetics have been published to date: one uses microsatellite markers to infer population structure (Davis et al., 2008) and two are based on the hyper-variable region (D-Loop) of the mitochondrial control region to infer genetic diversity and historical demography (Hernández-Ardila et al., 2021;Bender et al., 2023).In an era of constantly advancing sequencing technologies, the next logical step in leopard seal research is to apply genomic approaches to complement and enhance existing genetic data sets.Therefore, we present a high-quality reference genome for the leopard seal and evaluate its quality in comparison to the other four published genomes of the Monachinae subfamily.Organic tissue was collected during a necropsy of a stranded adult male leopard seal (Leopard Seal ID #: HLNZ-N013). The individual was found on August 5th, 2023, at Kaitorete Beach, Christchurch, New Zealand (43° 50' 0.781" S, 172° 34' 43.0" E) and was not yet decomposed at the time of necropsy (August 6th, 2023). Tissue samples were stored in RNA later and frozen at -20°C. Sample collection was conducted under the permit number 63499-MAR (New Zealand), and samples were transported to the United States under permit NMFS permit #26767 (USA, valid from 07/11/2023 to 06/30/2028).Genomic DNA was extracted from kidney tissue using the Qiagen Genomic DNA Extraction Kit, following the manufacturer's protocol. DNA quantity and quality were assessed using a Genomic DNA Screen Tape of a TapeStation System (Agilent) to ensure the presence of high-molecularweight DNA, required for long-read HiFi sequencing (Supplementary Figure S1). Library preparation was performed using the SMRTbell Express Template Prep Kit 2.0, and sequencing was conducted on a PacBio Sequel IIe system with five SMRT cells. The raw reads were assembled using Hifiasm v0.15.4-r347 (Cheng et al., 2021). The resulting contigs were queried against the NCBI nucleotide (nt) database using BLAST to identify potential contaminants (Camacho et al., 2009). These results were analyzed with BlobTools v1.1.1 (Laetsch and Blaxter, 2017), and contaminants were identified and removed from the assembly. To further refine the assembly, haplotigs, and redundant contig overlaps were removed using purge_dups v1.2.5 (Guan et al., 2020). All laboratory procedures, including DNA extraction, HiFi sequencing, genome assembly, and contaminant analysis (including graphical visualization with BlobTools), were carried out by Cantata Bio/Dovetail Genomics (CA, USA).Cantata Bio/Dovetail Genomics conducted a quality assessment of the genome assembly. This included the calculation of N50, L50, N90, and L90 statistics, and a genome completeness analysis using BUSCO version 4.0.5 and the eukaryota_odb10 database (2020-09-10, number of species: 70, number of BUSCOs: 255; Manni et al., 2021). We calculated additional statistics such as GC content and the length of the longest contig (in bp) using the stats.sh script of bbtools v39.06 (Bushnell et al., 2017). To receive an additional level of genome completeness, we also repeated the BUSCO analysis with the carnivora_odb10 lineage dataset (Creation date: 2024-01-08, number of genomes: 12, number of BUSCOs: 14502; BUSCO version 5.7.0; Manni et al., 2021). Quality statistics of our H. leptonyx assembly were compared to published reference genomes of the Monachinae subfamily (L. weddellii: GCA_000349705.1, M. angustirostris: GCA_029215605.1, M. leonina: GCA_011800145.1, and N. schauinslandi: GCA_002201575.2). We conducted further quality and completeness assessment via synteny analyses with the high-quality reference genomes of closely related species. Our reference genome was aligned to the chromosome-level genome of a N. schauinslandi male using minimap2 v.2.28-r1209 (Li et al., 2018). Although minimap2 is generally used as a sequence mapping tool, the parameters were adjusted for pairwise whole genome alignment as suggested in the current software manual. We also aligned our reference genome to the scaffold-level genome of a female individual of the more closely related species M. angustirostris (Accession #: GCF029215605.1). Synteny plots were created using the JupiterPlot pipeline (https://github.com/JustinChu/JupiterPlot; 2024), which implements the plotting tool Circos (Krzywinski et al., 2009). For both plots, we chose to display a maximum number of 50 scaffolds/contigs (maxScaff) while showing only contigs that are larger than 0.05% of the total length of the leopard seal reference genome. The minimum scaffold/contig size of the reference species (m) was set to 5,000,000 bp, ensuring the inclusion of all chromosomes of N. schauinslandi and reducing the number of scaffolds for M. angustirostris to 17. Maximum gap length (maxGap), minimum mapping quality (MAPQ), and minimum bundle size (minBundleSize) were set to 20,000 bp, 55, and 100,000 bp, respectively. Finally, we assessed the amount of sequence repeats and identified repeat families in the reference genome using the software packages RepeatModeler v.2.0.5 and RepeatMasker v.4.1.2 with de novo default settings (Smit, Hubley and Green RepeatMasker Open-4.0. 2013-2015 http://www.repeatmasker.org).The leopard seal reference genome consists of 203 contigs with a total size of 2.4 Gbp. The genome assembly shows an average coverage of 86X and the number of contigs with a length sum corresponding to half of the genome (L50) being 9. The N50 value, which represents the sequence length of the shortest contig representing 50% of the genome, is 99.45 Mbp. The assembly has a GC content of 41.6%, and the BUSCO completeness scores are 94.9% and 98.2% using the eukaryote and carnivora datasets, respectively (Supplementary File S1 and File S2). Based on these quality criteria, it is comparable to other pinniped reference genomes (Table 1). Moreover, it is an improvement over the only other currently available Lobodontini genome (L. weddellii; Noh et al., 2022) in terms of L50 and N50 values, and overall BUSCO completeness (Table 1). Thus, this leopard seal reference genome is likely the most complete Lobodontini genome to date. Approximately 35.91% of the reference genome is composed of repetitive elements. While some contigs consist almost exclusively of repetitive sequences, the proportion of repeats across the nine largest contigs ranges between 33.06% and 35.04% (L50 = 9; Supplementary Figure S2). These values are similar to the proportions reported for other phocid genomes (e.g. M. leonina = 41.51% Kim et al., 2020;Phoca largha = 35.83% Park et al., 2018). Out of the 35.91% repetitive sequences, 9.79% are retroelements (LINEs = 9.23%; LTRs = 0.55%), 0.43% are DNA transposons, 1.16% are simple repeats, 0.22% elements have low complexity and 24.32% are unclassified (Supplementary File S3). The high quality and completeness of this reference genome suggests that it will be eminently suitable for a variety of applications, both within and among species.The genome assembly data and raw reads are deposited at the Genebank repository of the NCBI database (http://www.ncbi.nlm.nih.gov/) under the BioProject and BioSample ID PRJNA1194539and SAMN45188358, respectively. The obtained genome assembly was submitted and registered under the NCBI GenBank accession number JBJQNM000000000. Raw long reads are publicly available at the NCBI Short Read Archive (SRA) under accession number SRR31619110.The whole genome alignment reveals that the H. leptonyx genome assembly covers 94.85% and 93.86% of the N. schauinslandi and M. angustirostris reference genomes, respectively. In the synteny analysis, 45 contigs of the H. leptonyx reference genome are assigned to the 18 N.schauinslandi chromosomes, of which three contigs map to the Y chromosome (scaffold IDs: ptg0036, ptg0059, and ptg0142). With the same parameters, 42 of these 45 contigs also match the 17 longest scaffolds of the M. angustirostris reference genome (Figure 1). The difference in the total number of contigs revealing homology is due to the three contigs mapping to the Y chromosome of N. schauinslandi, as the M. angustirostris genome assembly belongs to a female.Furthermore, either individual long contigs or merged contigs correspond to each chromosome of N. schauinslandi (Figure 1). This is another indication of the high level of genome completeness, as chromosome numbers among pinnipeds are highly conserved (2n = 34 to 2n = 36, Arnason, 1974;Beklemisheva et al., 2020). Long-term observations of this apex predator are largely limited by the remoteness of the leopard seal's habitat (i.e. pack-ice surrounding Antarctica) and their solitary nature. The resulting knowledge gaps about their basic biology and ecology prevent the effective incorporation of the species into ecosystem models. By presenting the first high-quality reference genome of the leopard seal, this study makes a significant step toward closing these gaps, providing a basic tool for future genomic analyses of the species' molecular ecology and evolutionary history, as demonstrated for other pinnipeds (Yakupova et al., 2022;Hoffman et al., 2024;Hauser et al., 2024). In particular, genomic resources enable analyses of the leopard seal's fine-scale population structure and kinships and can be used for a more robust assessment of its past and present effective population size. Furthermore, this leopard seal reference genome serves as a valuable tool for examining genetic adaptations and the evolution of key traits. It also facilitates investigations on its adaptability and potential to respond to rapid environmental change, especially as leopard seals appear to be inherently vulnerable to ice loss in many areas of their global distribution (Bender et al., 2023;Borras-Chavez et al., 2024), and can be negatively affected by anthropogenic actions such as overharvesting at the lower levels of the food chain (Forcada et al., 2009).