Background & Summary

The finless porpoise (Neophocaena spp.) is a group of small-sized, toothed whales that are mainly distributed in southern and eastern Asia. Their distribution includes the coastal waters of the western Pacific Ocean, Indian Ocean, Sea of Japan, and they also appear in the Bohai Sea, Yellow Sea, East China Sea, South China Sea, and middle and lower reaches of the Yangtze River in Chinese waters1,2. Since Cuvier first named the species Delphinus phocaenoides in 1829, the taxonomy and nomenclature of the finless porpoise have been controversial3,4. For decades, the finless porpoise was considered to be a single species consisting of three subspecies5,6,7, until Wang and Jefferson et al. concluded that the genus Neophocaena can be divided into two separate species based on their morphological and genetic characteristics, including the Indo-Pacific finless porpoise (N. phocaenoides) and the narrow-ridged finless porpoise (N. asiaeorientalis). The narrow-ridged finless porpoise can also be divided into two subspecies that include the Yangtze finless porpoise (N. a. asiaeorientalis) and the East Asian finless porpoise (N. a. sunameri)8,9, and this classification has been generally accepted. In 2018, Zhou et al. performed de novo genome sequencing of the Yangtze finless porpoise and re-sequenced three geographic populations in Chinese waters to investigate the freshwater adaptation mechanisms of the Yangtze finless porpoise10. Their results found that the genetic differentiation between the Yangtze finless porpoise and East Asian finless porpoise reached interspecific level, which supports their classification as independent species10.

With conservation, The IUCN Red List of Threatened Species categorized the Yangtze finless porpoise as “critically endangered” in 201311, and the narrow-ridged finless porpoise as “endangered” in 201712. However, the East Asian finless porpoise was not listed separately. The East Asian finless porpoise was listed in the Second Class of the National Key Protected Wild Animals List in China announced on February 5, 2021. Similar to other small cetaceans throughout the world, the East Asian finless porpoise population faces many critical factors, such as marine environment pollution, fishing injury, loss of important habitat, and decline of fish resources under the dual influence of global climate change and human activities13. Ultimately, the prospect of East Asian finless porpoise increasing in population is not optimistic, and it is extremely urgent to explore more conservation efforts for this species.

With regards to field research, East Asian finless porpoises are mainly found in the temperate waters of the west coast of the Pacific Ocean. For example, the coastal waters from the Taiwan Strait through the East China Sea north to the Yellow Sea and the Bohai Sea in China, as well as the waters of Korea and Japan7 (Fig. 1). Their wide distribution area makes it incredibly difficult to conduct systematic survey assessments for the entire population. Therefore, their total population size has not been accurately reported, which is unfavorable for their conservation. Because of these issues, population numbers are only recorded in local areas. Yoshida, Amano and Shirakihara assessed the population size of East Asian finless porpoise in the Ariake Sound/Tachibana Bay, Chiba/Sendai Bay and the Inland Sea of Japan, and estimated the population size to be 7,572 in the Inland Sea of Japan14,15,16. Population evaluation of the East Asian finless porpoise in China is only found in the Bohai Sea17.

Fig. 1
figure 1

Location distribution and photograph of the East Asian finless porpoise, N. a. sunameri. (a) A natural distribution map of N. a. sunameri and the sampling site (red star) of the study. (b) N. a. sunameri photographed in Penglai Polar Ocean World, Penglai, Shandong, China.

Full size image

Compared to Yangtze finless porpoise, very few systematic studies have explored the molecular biology, ecology, acoustics and feeding behavior of the East Asian finless porpoise over the past 40 years. Ruan et al. sequenced and compared the renal transcriptomes between the Yangtze finless porpoise and the East Asian finless porpoise to investigate the mechanism of osmotic pressure regulation with adaptation to their different habitats18. Additionally, Li et al. used a single hydrophone to record and analyze the echolocation signals of East Asian finless porpoises in Liaodong Bay and conducted a comparative study to Yangtze finless porpoises19. Further, Dong et al. concluded that the migration pattern of the East Asian finless porpoise population is mainly related to the migratory distribution of its preferred fish20, and finless porpoises have a broad diet that largely consists of fish, shrimp, and cephalopods21,22. Although these studies help understand finless porpoise migration behaviors and adaptation, more research needs to focus on improving conservation efforts for the East Asian finless porpoise.

In China, conservation research on the East Asian finless porpoise is less enthusiastic than that on the Yangtze finless porpoise. There is a serious lack of basic research on the East Asian finless porpoise, especially on its current population size, distribution characteristics, migration patterns, and key habitats. To date, the population size and distribution pattern of the East Asian finless porpoise in China are still not systematically known, and little is known about its key habitats. Consequently, its conservation biology research should receive more attention because it is an endangered marine mammal that is also listed as a second-class key protected wild animal in China.

The goal of this study was to assemble a gapless genome for the East Asian finless porpoise to aid in the conservation of this species. Here, we report a gapless cetacean genome that was generated through combining PacBio HiFi long reads and Hi-C sequencing data. We sequenced and analyzed the genome of the East Asian finless porpoise at the chromosomal level to gain a deeper understanding of its genetic background and evolutionary characteristics. The assembled genome size is approximately 2.50 Gb with a contig N50 of 84.69 Mb and scaffold N50 of 122.40 Mb. A total of 52 contigs were anchored onto 23 chromosomes (21 autosomes, X and Y chromosomes), and one mitochondrion chromosome. This genome contained 7 gapless assemblies of chromosome 4 (147 Mb), chromosome 9 (108 Mb), chromosome 11 (102 Mb), chromosome 16 (86 Mb), chromosome 18 (80 Mb), chromosome 21 (36 Mb), and chromosome X (125 Mb). Consequently, only 28 gaps were retained for next step filling in our assembly results. As the telomere-to-telomere assembly of human genome published this year23, ultra-long (>100-kbp) nanopore reads can be enable to span complex repeats and complete assemblies of the centromeres and telomeres. Gene annotation yielded 22,814 protein-coding genes and 97.31% of the predicted genes were annotated in publicly available biological databases, including NR, GO, KOG, KEGG, TrEMBL, Interpro and Swissprot. This high-quality assembled genome will provide rich research resources for conservation biology and phylogenetic studies on the East Asian finless porpoise, as well as research on genetic differentiation and adaptive evolution of other small toothed whales, like the Yangtze finless porpoise.

Methods

Sample collection

A muscle sample was collected from a male specimen of East Asian finless porpoise that died in the Yellow Sea near Lianyungang City, Jiangsu Province, China, in 2019 (Fig. 1). No ethical issues were considered in this study. The muscle sample was washed 3 times with Phosphate buffer saline (PBS), quickly frozen in liquid nitrogen, and stored at −80 °C until DNA extraction.

WGS library construction and genome size estimation

DNA was extracted from muscle specimen of the East Asian finless porpoise using MZ 1.3 (hypervariable minisatellite probe), as well as locus-specific minisatellite probes (g3, MS1 and MS43). For the short insert WGS library, DNA was sheared into fragments between 50 to 800 bp using a Covaris E220 ultrasonicator (Covaris, Brighton, UK) according to the manufacturer’s instructions. Fragments between 300 to 400 bp were selected to generate a single-stranded circular DNA library. The DNA library was sequenced on a MGISEQ-2000 platform. A total of 232.16 Gb of raw short reads were generated and 182.87 Gb of clean data were retained after adaptor removing and low-quality reads filtering by SOAPnuke (v2.0)24 with parameters “-n 0.01 -l 20 -q 0.1 -i -Q 2 -G 2 -M 2 -A 0.5” (Supplementary Table S1).

We used KmerGenie (v1.7051)25 to estimate the genome size with varied k-mer sizes from 21 to 121 (Fig. 2a). According to the smooth curves of estimated genome size, we obtained the predicted optimal k-mer size of 91 and the predicted genome size of 2,475,638,739 bp (Fig. 2b). The predicted genome size of the East Asian finless porpoise is consistent with that of the Yangtze finless porpoise (2.49 Gb) found in a previous study10.

Fig. 2
figure 2

K-mers analysis to estimate the genome size of the East Asian finless porpoise. (a) 182 Gb of high-quality data was used to generate 16 different k-mers depth distribution curve frequencies by KmerGenie. The k-mers value was automatically set by the software from 21 to 121. The x-axis indicates k-mer size, while the y-axis is the number of genomics k-mers at that k-mer size. (b) 91-mer depth frequency distribution. The x-axis is depth (X), while the y-axis is the proportion that represents the frequency at that depth divided by the total frequency of all depths. The genome size was estimated using the following formula: Genome size = K-mer num/Peak depth. The peak depth is approximately 28 and the estimated genome size is 2,475,638,739 bp.

Full size image

PacBio library preparation, sequencing, and de novo assembly using HiFi reads

DNA was extracted from the same muscle specimen using a QIAGEN Blood & Cell Culture DNA Midi Kit following the manufacturer’s instruction (QIAGEN, Germany). After DNA extraction, two sequencing libraries were prepared according to the “Using SMRTbell Express Template Prep Kit 2.0 With Low DNA Input” protocol from PacBio with an insert size of approximately 20 kb (Pacific Biosciences, USA). The libraries were then sequenced on a PacBio Sequel II SMRT cells in circular consensus sequence (CCS) mode. A total of 5 SMRT Cells were sequenced. 2,397 Gb subreads were processed using the CCS algorithm of SMRTLink (v8.0.0)26 with parameters “–minPasses 3 –minPredictedAccuracy 0.99 –minLength 500”, yielding 154 Gb of PacBio’s long high-fidelity (HiFi) reads (Supplementary Table S1). With the HiFi reads, the primary contigs were assembled by Hifiasm (v0.15.1)27 with default parameters. After, the Purge-Haplotigs28 program was used to remove redundant sequences with parameters “-j 80 -s 80 -a 30”, which yielded a contig assembly with a size of approximately 2.50 Gb and contig N50 of 84.69 Mb (Table 1).

Table 1 Statistics of assembly.
Full size table

Hi-C library preparation, sequencing, and chromosome anchoring

The same muscle specimen was fixed with 1% formaldehyde for 10–30 min at room temperature to coagulate proteins that are involved in chromatin interaction in the genome. The restriction enzyme Mbo I (NEB, Ipswich, USA) was then added to digest the DNA, and fragments with flat or sticky ends were obtained. The ends were flattened and repaired, and then labeled with biotin. The inter-match fragments were ligated with T4 DNA ligase (Thermo Scientific, USA) to form a loop. Proteins that connected the DNA fragments were then digested to obtain the crosslinked fragments, and the clip was interrupted again using ultrasound. A Hi-C library was made by capturing the biotin with magnetic beads and sequenced on a MGISEQ-2000 instrument. A total of 219.2 Gb of clean data were obtained from 263.87 Gb of sequencing data using software SOAPnuke (v2.0)24 with parameters “-n 0.01 -l 20 -q 0.1 -i -Q 2 -G 2 -M 2 -A 0.5” (Supplementary Table S1).

To anchor contigs onto chromosomes, the Hi-C clean data were mapped to the assembled contigs using BWA (v0.7.12)29, and then erroneous mappings (MAPQ = 0) and duplicates were filtered by the juicer pipeline (v1.5)30 to obtain the interaction matrix. Following, approximately 625.70 Mb reads (~77.77%) were used to anchor the contigs into chromosomes with 3D-DNA pipeline (v180,922)31. And 3D-DNA pipeline31 was used to remove select short contigs using default parameters. The Hi-C contact maps were then reviewed with JUICEBOX Assembly Tools (v2.15.07)30. These processes generated a final genome assembly, where the genome size was approximately 2.50 Gb and contig N50 was 84.69 Mb. Remarkably, 52 contigs were linked onto the 21 autosomes, two sex chromosomes, and one mitochondria sequence (Fig. 3, Tables 1 and 2).

Fig. 3
figure 3

Genome-wide Hi-C heat maps of the East Asian finless porpoise genome.

Full size image
Table 2 Chromosomes Length of East Asian finless porpoise.
Full size table

Identification of Y chromosome sequences

Generally, sequence assembly on the Y chromosome is a challenge due to its complex repetitive nature. Here, we assembled all the PacBio HiFi reads into contigs by Hifiasm27 software, and then filtered the redundant sequences using Purge-Haplotigs28 software. Finally, we anchored the non- redundant contigs into scaffold with Hi-C data. In order to identify the Y sequence, we mapped the scaffold sequences onto the Y chromosome of Tursiops truncates32. Additionally, we mapped the contig sequences of the East Asian finless porpoise genome into Y chromosome of Tursiops truncates genome using Ragtag (v2.1.0)33 with default parameters. Based on the above two methods, we obtained a candidate Y sequence with a high degree of similarity. We selected the Y sequence assembled from the first method. The newly assembled Y chromosome sequence is 11.02 Mb in length and contains 82 intact protein-coding gene models. Of these 82 genes, seven genes were identified as TSPY genes. In humans, these four genes (SRY, TSPY1, TSPY3 and ZFY) are linked to the Y chromosome34,35,36,37,38,39. Accordingly, this assembled Y chromosome sequence is highly reliable.

Repeat annotation

Two strategies including de novo and homolog methods were used to annotate repeat elements. De novo repeats were identified by RepeatModeler (v1.0.4)40 and long terminal repeats were annotated by LTR-FINDER (v1.0.7)41. DNA and protein transposable elements (TEs) were detected by RepeatMasker (v4.0.7)42 and RepeatProteinMasker (v4.0.7), respectively, based on Repbase database43. Tandem repeats were performed by Tandem Repeat Finder (v4.10.0)44. We obtained approximately 1.05 Gb (~42.23%) of repetitive sequences (Supplementary Table S2), which were similar to the Yangtze finless porpoise10, and 38.63% belonged to LINE subfamily (Supplementary Table S3).

Protein-coding genes prediction and functional annotation

To predict genes, we generated 32 RNA-seq samples from blood tissue of the Yangtze finless porpoise specimen (Supplementary Table S4)45,46. These reads were then aligned to East Asian finless porpoise genome using Hisat2 (v2.1.0)47 with the following parameters: –sensitive –no-discordant –no-mixed -I 1 -X 1000 –max-intronlen 1000000. The aligned reads were assembled using Stringtie (v1.3.5)48 with the following parameters: -f 0.3 -j 3 -c 5 -g 100 -s 10000. Subsequently, TransDecoder (v5.5.0) (https://github.com/TransDecoder/TransDecoder) was used to identify the coding sequence with default parameters. Gene models for the East Asian finless porpoise were also predicted by Augustus (v3.2.1)49 for de novo annotation. Homologous proteins of eight reference species were downloaded from common databases. Data for the Bowhead whale50 was downloaded from the Bowhead Whale Genome Resource (http://www.bowhead-whale.org). Common minke whale51, Beluga whale52, Yangtze River dolphin53, Yangtze finless porpoise10, Killer whale54, Bottlenose dolphin54 and Sperm whale were downloaded from the National Center for Biotechnology Information (NCBI). GeMoMa (v1.8)55 was used to search coding structures based on transcriptome data and homologous proteins. A total of 22,814 coding genes (36,167 transcripts) were predicted (Supplementary Table S5). All protein coding genes were supported by at least one prediction method (Supplementary Table S6). The final gene set was functionally annotated by mapping against KEGG56, Swiss-Prot57, TrEMBL57, KOG58, InterPro59 and NR (NCBI Non-redundant protein) databases using BLAST (v2.2.26)60 with an E-value threshold of 1E-5. The protein domains and motifs were annotated using InterProScan61. GO Ontology (GO)62 was obtained from the InterProScan61 results in this study, and 97.31% of the 22,814 proteins were annotated by at least one database (Supplementary Table S7). Of these functional proteins, 16,455 (~72%) were supported by all five databases (Fig. 4).

Fig. 4
figure 4

Gene function annotation results in a statistics Venn diagram using five public databases: NR, InterPro, KEGG, SwissProt and KOG.

Full size image

Data Records

The DNA sequence reads of East Asian finless porpoise (Experiment of DNA sequencing data from genome survey library: SRR2104715463; Experiments of DNA sequencing data from Hi-C library: SRR2076093564-SRR2076093665; Experiments of DNA sequencing data from PacBio HiFi library: SRR20997931-SRR2099793566,67,68,69,70) have been deposited in the Sequence Read Archive (SRA) under project number SRP38952971. The Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession JANJGR00000000072. Files of the assembled genome, gene structure annotation, repeat predictions and gene functional annotation of East Asian finless porpoise were deposited at Figshare database under DOI code73.

Technical Validation

Evaluation of the genome assembly

By comparing the assembled metrics of the East Asian finless porpoise to the other cetacean species, our assembly substantially improved because of increased contig and scaffold lengths, which indicates that our assembly is highly contiguous. Our gapless genome assembly increased the contiguity metrics 941-fold (by contig N50) or 921-fold (by the number of contigs) compared to a previously reported Yangtze finless porpoise assembly10. Among the public cetacean genomes, our assembly had the longest contig N50 length and smallest gap number, which suggests that our East Asian finless porpoise genome is high quality (Table 3).

Table 3 Comparison of the East Asian finless porpoise genome with previously published cetacean genomes.
Full size table

To assess the completeness of our East Asian finless porpoise genome, we performed an analysis using Benchmarking Universal Single-Copy Orthologs (BUSCO, v5.1.0)74 with the mammalia_odb10 database. The results showed that 95.4% of the expected mammalian genes (including 93.7% single and 1.7% duplicated ones) had complete gene coverage, and 1.1% were identified as fragmented, respectively. However, 3.5% were considered missing in our East Asian finless porpoise genome. Still, the complete evaluation of the East Asian finless porpoise genome is more superior than other current public cetacean genomes (Table 3).

To evaluate the telomere sequences assembled in the East Asian finless porpoise genome, we used the Telomere Identification toolkit (Tidk, v0.2.0) (https://github.com/tolkit/telomeric-identifier) to search telomere sequences (TTAGGG) along with the genome sequence. From the results, 23 chromosomes detected at least one side of telomere sequences, such as Chr5 and Chr11. Individual sequences were identified with partial telomere sequences, which should be further investigated and optimized.

To compare the genome consistency between the East Asian finless porpoise and the Yangtze finless porpoise, we used MuMmer (v4.0.0)75 to identify similar regions with parameters “–mum -c 500 -l 40” at the genome level. Additionally, we also used BLAST60 and WGDI (https://github.com/SunPengChuan/wgdi) software to search the synteny blocks with at least ten gene pairs at the gene level. These two analyses revealed that the two genomes are highly consistent (Fig. 5).

Fig. 5
figure 5

Comparison of sequence synteny between the East Asian finless porpoise and the Yangtze finless porpoise. (a) MUMmer was used to identify similar regions between the East Asian finless porpoise and Yangtze finless porpoise genome sequences. (b) WGDI was used to detect syntenic blocks between the East Asian finless porpoise and Yangtze finless porpoise gene pairs. The x-axis is the chromosome scale of East Asian finless porpoise genome, and the y-axis is the contig scale of Yangtze finless porpoise genome.

Full size image

Evaluation of the gene annotation

We performed BUSCO74 analysis with the mammalia_odb10 database to assess the completeness of the coding sequences for the East Asian finless porpoise. The results showed that 97.9% of the expected mammalian genes (including 96.7% single and 1.2% duplicated ones) had complete gene coverage, and only 0.5% were identified as fragmented, respectively. However, 1.6% were considered missing in our East Asian finless porpoise genome. Compared to other complete evaluations of protein-coding genes, our East Asian finless porpoise has a high degree of integrity (Table 4).

Table 4 Comparison of the East Asian finless porpoise genes with other representative cetacean gene sets.
Full size table