1 Introduction

At least 20000 species of fish, one of the largest groups of vertebrate animals, colonize virtually in freshwater environments as well as tropical, temperate, and polar saltwater environments. Olfaction is critical to many animals, including fish, to facilitate migration, avoidance of predators, reproduction, and feeding (Laberge and Hara, 2001). For example, fish inhabiting in rocky reefs tend to be weak swimmers, but their strong sense of smell helps them to avoid predators (Wei et al., 2014). The recent availability of taxonomically diverse fish genomic sequences makes it possible to conduct evolutionary analyses of olfactory receptor (OR) repertoires across fish species. The majority of studies on the genetic basis of fish OR repertoires have focused on the identification and characterization of individual genes, leaving the evolutionary history and the roles of fish OR genes in ecological adaptation relatively unexplored.

The order Lobotiformes (class Actinopterygii) contains three families (Datnioididae, Hapalogenyidae, and Lobotidae) and several economically important species (Betancur et al., 2013; Betancur et al., 2017). Lobotiformes is primarily composed of rocky reef species with limited or no long-distance migration. Hapalogenys analis belongs to the family of Hapalogenyidae and is a typical rocky reef fish native to East Asia, inhabiting in the coastal waters of China, Japan, the Korean Peninsula, and some countries in Southeast Asia. H. analis has high economic values, with delicate meat containing a high content of amino acids, and is used in marine ranching and offshore cage aquaculture in Korea, Japan, and China (Xu et al., 2010). As a predator that is active at night, what factors help it feeding and odorant-oriented movement? Is it the olfactory system? Genomic analyses can provide relevant information to these questions. However, the only high-quality genomic sequence available for a member of Lobotiformes is of Datnioides undecimradiatus (GenBank assembly accession: GCA_008933995.1). This lack of relevant genomic resources hinders the continued study of the taxonomy, evolution, and biology of H. analis as well as other Lobotiformes.

In this research, a chromosome-level genome assembly of H. analis was studied utilizing Hi-C and PacBio longread sequencing data. The genome was assembled to provide insight into the evolutionary history of the gene families related to OR repertoires of H. analis and the mechanism by which this species adapts to rocky reef environments. Through comparative genomics, we identified the genetic basis of enhanced olfaction function and weak swimming capability in rocky reef fish. This genomic assembly will be crucial for studies on the ecology, biology, and evolutionary history of H. analis, as well as other species of Lobotiformes.

2 Materials and Methods 2.1 Sample Collection and Extraction of DNA

In September 2020, an adult female fish weighing approximately 650 g was sampled from Zhoushan (Zhejiang, China) fishery located in southeast of the Yangtze Estuary. DNA degradation was prevented by storing the collected muscle tissue in liquid nitrogen prior to DNA extraction. Genomic DNA was isolated from muscle tissues using Blood & Cell Culture DNA Mini Kit (QIAGEN, Cat. No. 13343), which is based on the phenol/chloroform method (Sambrook et al., 1989). Total RNA was extracted using TRIzol reagent (Invitrogen). The quality and concentration of extracted DNA/RNA were assessed using a Pultton DNA/Protein Analyzer (Plextech), and their integrity was further evaluated on 1% agarose gel stained with ethidium bromide (EB). The extracted DNA/RNA were stored at −80℃ for further analyses. All tissue collection and DNA/RNA extraction were conducted following the ethical regulations provided by the Institutional Animals Care and Use Committee of Zhejiang Ocean University. Only high-quality DNA was used for PacBio Sequel (for genome assembly), Illumina (for genome survey and genome correction), and Hi-C (for assisted genome assembly) sequencing. Other tissues, including the eyes, liver, and spleen, were utilized for RNA extraction. Fresh tissue was used for Hi-C sequencing and the rest of the samples were rapidly frozen in liquid nitrogen for 1 h and then stored in an ultra-low temperature refrigerator at −80℃. The resultant sequencing data were used for genome annotation.

2.2 Construction of the Library and Sequencing of the Genome

Genomic sequencing was carried out using the Illumina NovaSeq-6000 platform (Illumina, United States), and longread generation was carried out using the PacBio Sequel Ⅱ platform (Pacific Biosciences). The sequencing library generated by the Illumina platform were used for genome size estimation, genome assembly correction, and genome assembly evaluation. According to Illumina's standard procedure, a pair-end library was created utilizing a 300 bp insert size. The reads with more than 10% N bases or lowquality bases ≤ 5, adapter sequences, and duplicated sequences were removed. For the PacBio platform of long reads, a SMRTbell (Single molecule real-time, SMRT) library was established with a fragment size of 20 kb using the SMRTBell template preparation kit 1.0 (PacBio) following the manufacturer's instructions. The library was sequenced with the PacBio Sequel Ⅱ system in Circular Consensus Sequence (CCS) mode, and data from one SMRT cell were generated.

2.3 Assembly of the Genome and Estimation of Genome Size

Illumina short reads were subjected to K-mer analysis in order to estimate the repeat content, heterozygosity, and size of the H. analis genome (Liu et al., 2013). Genome size was calculated as follows:

$ G{\text{(enome)}} = \frac{{{\text{Number of }}K{\text{ - mers}}}}{{{\text{Average depth of }}K{\text{ - mers}}}}, $

where

$ {\rm{Number\;of\;}} K-{\rm{mers}} =\\ \;\;\;\;\;\; {\rm{Total}}\;K-{\rm{mer\;count–Abnormal\;}} K-{\rm{mer\;count}}, $

and depth of K-mers is peak depth. Canu was used to assemble PacBio long-reads, which were employed to polish the assembly data following removal of redundant data and haplotigs (Koren et al., 2017).

A two-step quality control process was adopted because of the errors intrinsic to long-reads. First, long-reads were used to correct the assembly with the GenomicConsensus (version 2.3.3) Arrow package, with minCoverage set to 15 (Chin et al., 2013). Second, the assembly was subjected to two rounds of quality control utilizing Illumina shortreads with PILON (version 1.2.3) (Walker et al., 2014). Following quality control, Purge_haplotigs was utilized to remove redundancy, according to sequence similarity and read depth (Roach et al., 2018).

2.4 Chromosome-Level Genome Assembly with Hi-C Library

The method of Kajitani et al. (2014) was used to assemble a chromosome-level genome utilizing a Hi-C library. An Illumina Novaseq 6000 platform was used to sequence the library, and bowtie was used to map the reads to the quality-controlled genome, using standard settings (Langmead, 2010). Both read ends were aligned independently, and only read pairs exhibiting unique alignment at each end were selected for further analysis. Lachesis (version 1.03) was used to assemble the chromosome-level genome with corrected contigs and Hi-C reads, using standard settings (Burton et al., 2013). Juicebox (version 1.5) was used to manually adjust the positions and directions of small contigs, according to the degree of contig interaction (Durand et al., 2016).

Two methods were employed to evaluate the consistency and accuracy of the assembled chromosome-level genome. First, Blasr (version 5.3.3) was used to map the longreads to the genome (Chaisson and Tesler, 2012), and BWA-MEM (version 0.7.10) was used to map the short-reads (Li et al., 2010), using default settings. Second, genomic assembly completeness was estimated using Benchmarking Universal Single Copy Orthologues (BUSCO version 2.0), based on single copy homologs (Simão et al., 2015).

2.5 Gene Model Prediction and Functional and Repeat Annotation

The repeat elements were annotated by homology-based and de novo prediction methods. The homolog prediction was conducted with Repbase database (https://www.girinst.org/repbase/) (Jurka et al., 2005). Repeat Protein Mask (version 4.1.0) and RepeatMasker (version 4.0.7) were utilized to extract repeat regions (Tarailo-Graovac and Chen, 2009). Tandem repeats were extracted using Tandem Repeat Finder (version 4.09) (Benson, 1999) by ab initio prediction. LTR_FINDER (version 1.0.2) (Xu and Wang, 2007), Repeat Modeler (version 2.0) (Flynn et al., 2020), and Repeat Scout (version 1.0.6) (Price et al., 2005) were used for de novo identification of transposable elements (TEs). All TEs and repeats were combined into a repeat library to further verify the annotation results of repeat sequences.

Gene structure prediction used a combination of de novo prediction, homology-based prediction, and transcriptomebased strategy. Both GENSCAN (Burge and Karlin, 1997) and AUGUSTUS (version 2.7) (Stanke et al., 2006) were used for de novo predictions. For homology-based annotation, the complete protein sequences of five related species were downloaded from NCBI (the National Center for Biotechnology Information), including Acanthochromis polyacanthus, Anabas testudineus, Astyanax mexicanus, Amphiprion ocellaris, and Astatotilapia calliptera. TBLASTN (E value ≤ 1e − 5) was used to compare the genome sequence of H. analis with them. GENEWISE (version 2.4.0) was used to predict gene model according to the corresponding gene region on each alignment (Birney et al., 2004), and the intron/exon structure was defined (Doerks et al., 2002). Thirdly, transcriptome data were aligned to the assembled genome sequences using TopHat (version 2.0.11) to identify exon regions and splice positions (Kim et al., 2013). The alignment results were then used as input for Cufflinks (version 2.2.1) for genome-based transcript assembly (Ghosh and Chan, 2016). The gene sets predicted by various methods were integrated into a complete gene set using HiCESAP (Wuhan Gooalgene Co., Ltd., https://www.gooalgene.com/) and Maker (version 2.31.10) (Cantarel et al., 2008). Functional annotation of predicted genes was carried out using Blastx/n (E value ≤ 1e − 5), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000), TrEMBL (Boeckmann et al., 2003), InterPro, SwissProt, and the NCBI non-redundant database (NR). Blast2GO (version 5.2.5) was utilized to perform Gene Ontology (GO) annotation (Conesa et al., 2005).

Four kinds of non-coding RNA were annotated, including miRNA, tRNA, rRNA and snRNA. Micro-RNAs and ribosomal RNAs were identified by searching the Rfam database using the Infernal (version 1.1) (Griffiths-Jones et al., 2005; Nawrocki and Eddy, 2013). TRNASCAN-SE (version 1.3.1) was used to predict transfer RNAs (Lowe and Eddy, 1997), and rRNAs were predicted using BLASTN.

2.6 Selection Analysis and Genome Evolution Analysis

OrthoMCL (verison v2.0.9) was used to detect orthologous groups by aligning the protein sequences of Liparis tanakae, Anarrhichthys ocellatus, Pseudochaenichthys georgianus, Trematomus bernacchii, Notothenia coriiceps, Etheostoma spectabile, Sebastes schlegelii, Sebastes umbrosus, Epinephelus lanceolatus, Larimichthys crocea, Argyrosomus japonicus, Collichthys lucidus, Danio rerio, and Lepisosteus oculatus. MUSCLE (version 3.8.31) was used to multiply-align the single-copy orthologs shared between all 15 species, using default parameters (Edgar, 2004). Raxml (version 8.2.12) was used to perform multiple-sequence alignment and phylogenetic tree construction (Stamatakis, 2014). Divergence time was evaluated with r8s (version 1.7) and the MCMC Tree package of PAML (Yang, 2007). Calibration was carried out using the time of divergence between L. oculatus and D. rerio (245.1–335.6 Mya), which was acquired from TimeTree (http://www.timetree.org/) (Kumar et al., 2017). Gene contraction and expansion were evaluated using CAFÉ (version 4.0) (De Bie et al., 2006). The contracted and expanded gene families in H. analis, and the contracted gene families shared between five sedentary species (H. analis, A. ocellatus, S. schlegelii, S. umbrosus, and E. lanceolatus) were found to be enriched in several GO terms and KEGG pathways, and the Benjamini and Hochberg FDR correction was applied. Using the KEGG and GO annotations of H. analis as background, the significantly overrepresented KEGG pathways and GO terms were identified (corrected P value ≤ 0.05). H. analis was used as foreground branch, and other species for which phylogenetic trees were constructed were used as background branches. Positively selected genes (PSGs) in the H. analis genome were detected using single-copy orthologous genes. The branch-site model and the likelihood ratio comparison test in the CodeML program of the PAML software were used to test whether individual codon sites in the gene sequences were selective (determined by calculating the site dN/dS values, i.e., non-synonymous substitution rate/synonymous substitution rate), and the possible positively selected genes obtained were subjected to GO and KEGG enrichment analyses.

3 Results and Discussion 3.1 Genome Assembly and Assessment

Illumina Novaseq sequencing was conducted using a 300 bp library, resulting in the generation of 83.02 Gb of sequence data. According to the K-mer depth and number, we inferred that the genome was approximately 550 Mb in size. The genome size was then corrected to 543 Mb through removal of incorrect K-mers resulting from sequencing errors. The heterozygosity rate was 0.45% and the repeat rate was 26.78%. Approximately 120.67 Gb (about 222×) of PacBio sequencing data was obtained utilizing the Pac-Bio Sequel Ⅱ platform. The final genomic sequences were 544.26 Mb in length, consistent with the K-mer-based estimate. The genome contained 755 contigs, with contig N50 and scaffold N50 values of 3.43 Mb and 23.53 Mb, respectively. This genome was compared to the Actinopterygii database, and BUSCO analysis indicated that 97.27% (4459/4584; 94.20% for single-copy and complete BUSCO and 3.08% for duplicated and complete BUSCO) of the complete BUSCO genes were present in the assembled genome.

Hi-C scaffolding was used to anchor and orient the genomic contigs to each chromosome (Fig. 1). 60.235 Gb of clean data were generated by the Hi-C library, of which 54.34% was validly paired. Using Lachesis, we anchored 99.02% of the assembled sequences to 24 chromosomes, with lengths ranging between 10.36 Mb to 28.90 Mb (Fig. 2). The 24 chromosomes were clearly distinguished in the Hi-C heatmap, and the strong internal interactions indicated the high quality of anchoring. The final assembled genome was 539.00 Mb in length, with a scaffold N50 of 23.53 Mb and a contig N50 of 3.43 Mb (Table 1).

3.2 Gene Model Prediction and Functional and Repeat Annotation

A total of 132.74 Mb of repeat sequences were annotated, which accounted for 24.39% of the total genome. The predominant repeat elements were DNA transposons (12.63%), long terminal repeats (LTR, 5.40%), and long interspersed element (LINE, 4.79%). A combination of RNA-seq based method, homolog-based method, and de novo method was used to predict the presence of 21360 protein-coding genes. The average number of exons, exon length, intron length, CDS length, and gene length were 10.44, 237.66 bp, 1126 bp, 1704 bp, and 13105 bp, respectively. The predicted gene model statistics exhibited similar patterns across A. polyacanthus, A. testudineus, A. mexicanus, A. ocellaris, and A. calliptera (Fig. 3). The BUSCO evaluation of predicted genes indicated that 96.90% of the total orthologous genes were matched (93.10% of singlecopy and complete BUSCO, and 3.80% of duplicated and complete BUSCO). A total of 20787 genes, which accounted for 97.32% of all genes, were annotated within public databases. We also successfully annotated 511 miRNAs, 2904 tRNAs, 1897 rRNAs, and 753 snRNAs, with average lengths of 85, 75, 122, and 146 bp, respectively. The complete assembly and annotation results were visualized in a circus plot (Fig. 4).

3.3 Comparative Genomics

A comparative genomic assessment between H. analis and 14 other species (five sedentary fish and seven migrating fish) revealed that the H. analis genome contained 3584 contracted and 364 expanded gene families, compared to the most recent common ancestor of H. analis and Sciaenidae (L. crocea, C. lucidus) (Fig. 5). Further analysis indicated that gene families under expansion were significantly enriched in 4 KEGG pathways and 36 GO terms associated with olfaction, including olfactory receptor activity (corrected P-value = 8.43E − 11, GO: 0004984), G protein-coupled receptor activity (corrected P-value = 1.04E − 05, GO: 0004930), transmembrane signaling receptor activity (corrected P-value = 5.38E − 05, GO: 0004888), olfactory transduction (corrected P-value = 8.55E − 11, ko04740), and necroptosis (corrected P-value = 8.91E − 07, ko04217). Olfaction is crucial to an array of animal behaviors, including predator avoidance, mate selection, and feeding (Hara, 1975; Su et al., 2009; Bazáes et al., 2013; Hughes et al., 2018). Through regulation of the OR gene family, animals can detect and discriminate between a wide variety of odiferous compounds. In both fish and mammals, OR gene expression is the highest in olfactory epithelium (OE) sensory neurons, which are located within the nasal cavity (Vassar et al., 1993; Churcher et al., 2015; van der Linden et al., 2018; Cong et al., 2019). Odorant discrimination and, by extension, environmental adaptation are driven largely by the diversity of OR genes. Among placental mammals in particular, there exists considerable variation in the number of intact OR genes, with most species hosting a large number of OR pseudogenes (Niimura et al., 2014). This variation in the OR repertoire is primarily driven by frequent evolutionary losses and gains of genes through pseudogenization and duplication events (Vandewege et al., 2016; Hughes et al., 2018; Niimura et al., 2018).

The H. analis genome was found to contain contracted gene families related to motor function, including myosin complex (corrected P-value = 1.48E − 28, GO: 0016459), actin cytoskeleton (corrected P-value = 3.02E − 24, GO: 0015629), cytoskeletal motor activity (corrected P-value = 9.13E − 23, GO: 0003774), carbohydrate derivative binding GO terms (corrected P-value = 5.83E − 10, GO: 0097 367), vascular smooth muscle contraction (corrected P-value = 1.01E − 10, ko04270), cardiac muscle contraction (corrected P-value = 1.35E − 8, ko04260), thyroid hormone signaling pathways (corrected P-value = 1.93E − 7, ko04919), and arachidonic acid metabolism pathways (corrected P-value = 3.31E − 4, ko00590). Notably, the five sedentary species (A. ocellatus, E. lanceolatus, H. analis, S. schlegelii, and S. umbrosus) shared 19 contracted gene families associated with cellular aromatic compound metabolic process (corrected P-value = 0.0E + 00, GO: 0006725), nitrogen compound metabolic process (corrected P-value = 0.0E + 00, GO: 0006807), metabolic process (corrected P-value = 0.0E + 00, GO: 0008152), biosynthetic process (corrected P-value = 0.0E + 00, GO: 0009058), and primary metabolic process (corrected P-value = 0.0E + 00, GO: 0044238). These functions are critical for regulating energy metabolism, and there is a correlation between lower metabolic rates and lower DNA mutation rates. Specifically, the DNA mutation rate is related to environmental energy (Davies et al., 2004), metabolic rate (Martin and Palumbi, 1993), and life history characteristics (Bromham et al., 1996). Certain contracted gene families, for instance neutrophil extracellular trap formation (corrected P-value = 0.0E + 00, ko04613), may reflect the reduced energy requirements and locomotor activity of sedentary fish compared with migratory fish. Sun et al. (2011) found that the protein-coding mitochondrial genes of migratory fishes are under greater selective pressure than non-migratory fishes. In addition, several studies have shown that a lower metabolic rate may increase the tolerance of functional genes to mutations (Child, 1939; Zhu et al., 2020). For example, Cassidy et al. (2019) found that Drosophila can slow down development by reducing metabolic rate. This process inhibits gene expression disorder resulted from the deletion of miRNA-activating genes or the downregulation of miRNA expression, and allows the fly to develop a normal phenotype. This phenomenon can explain why sedentary fish develop normally despite accumulating more mutations than migratory fish.

A branch-site model identified a total of 277 genes undergoing positive selection (positive selection genes, PSGs). These PSGs were found to be significantly enriched in five KEGG pathways and 316 GO terms, including DNA binding (corrected P-value = 1.07E − 02, GO: 0003677), gene expression (corrected P-value = 1.82E − 02, GO: 0010467), response to stress (corrected P-value = 2.52E − 02, GO: 00 06950), DNA repair (corrected P-value = 2.52E − 02, GO: 0006281), and homologous recombination (corrected P-value = 8.20E − 03, ko: 03440). These functions are associated with macromolecule biosynthesis and environmental adaptation.

4 Conclusions

We utilized both Hi-C data and PacBio long reads to assemble the 539 Mb chromosome-level genome of H. analis, which represents the first reference genome of the family Hapalogenyidae. Owing to the intrinsically long Pac-Bio sequencing, the assembled genome had a scaffold N50 of 23.53 Mb and a contig N50 of 3.43 Mb. The assembled genome and genome annotation will be beneficial to the further study of Lobotiformes species, as well as to an array of breeding, conservation, and phylogenetic studies of the family Hapalogenyidae. Comparative genomic analysis revealed that OR gene families were expanded, while motor capacity and energy metabolism gene families were contracted in H. analis and other sedentary reef fishes. These results also help explain how H. analis adaptes to rocky reef environments.

Data Availability Statement

The raw sequencing data, including RNA-seq [SRA accession: SRR16077758], Illumina short reads [SRA accession: SRR16002978], Hi-C [SRA accession: SRR16002979], and PacBio data [SRA accession: SRR16002979], are available in the NCBI Sequence Read Archive (SRA) database, and can be accessed with the project ID PRJNA765047. Please contact the corresponding author for more annotation information.

Acknowledgements

This work was supported by the Province Key Research and Development Program of Zhejiang (No. 2021C02047), the Special Projects of Zhejiang Provincial Science and Technology Department (Nos. HYS-CZ-004, HYS-CZ-20 2208), and the 'San Nong Jiu Fang' Science and Technology Cooperation Project of Zhejiang Province (No. 2022 SN JF073).