The sample was captured from the Ross Sea (77◦050 S, 170◦300 E on CCAMLR Subarea 88.1), Antarctica, and then placed in a -80 ◦C freezer. The specimen was used for total RNA extraction and DNA isolation. The frozen muscle samples were anatomized, and the traditional phenol-chloroform method was used to produce high-quality and pure DNA for long-read sequencing on the PacBio Sequel platform (Pacific Biosciences, USA). The quality and quantity of genomic DNA were confirmed using a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and a Fragment Analyzer (Agilent Technologies, Santa Clara, CA, USA), respectively. Following the manufacturer’s guidance, total RNA was extracted from muscle and skin samples using a RNeasy Mini Kit (Qiagen). The quality and quantity of the RNA were confirmed using a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and a Fragment Analyzer (Agilent Technologies, Santa Clara, CA, USA).

The genomic DNA was used to construct the single-molecule real-time (SMRT) bell library and was sequenced using eight SMRT cells (Pacific Biosciences, Sequel™ SMRT Cell 1M v2) with a Sequel Sequencing Kit 2.1 (Pacific Biosciences, USA) and one 600-minute movie was captured for each SMRT cell.

The muscle tissue was used to construct the whole-genome sequencing (WGS) library via the TruSeq Nano DNA library kit (Illumina San Diego, CA, USA). We used 1 µg of DNA, which was sheared into 350–450 bp fragments using a Covaris E220. The end of each DNA fragment was repaired, and an ‘A’ base was added to the 3’ end to make a sticky end. After ligation, an Illumina adapter was ligated to both ends of the DNA fragments, and PCR amplification was performed on each sample. The constructed library was purified, and the final WGS library was sequenced on an Illumina NovaSeq 6000 platform (Illumina San Diego, CA, USA).

The Hi-C chromatin contact map was constructed using muscle tissues from the same individual for the high-quality chromosome-level assembly. The muscle was crosslinked with PBS and formaldehyde, and the chromatin was prepared. 800 ng of the chromatin were utilized to prepare the library after the chromatin sample had been normalized. The chromatin was collected and digested with a restriction enzyme. Using a Covaris system after the crosslink reversal, 200 ng of DNA was sheared. The sheared DNA fragments were end-repaired and ligated with an Illumina adaptor. The ligated DNA was verified using Streptavidin magnetic beads before being amplified via PCR to enrich the fragments. The insert size and the concentration of the constructed libraries were checked. Then they were sequenced using an Illumina NovaSeq 6000 platform (Illumina San Diego, CA, USA) with a 150 bp paired-end strategy.

Total RNA (1 µg) from the muscle and skin tissues was pooled for RNA sequencing. The pooled samples were sequenced using one SMRT cell v3 based on P6-C4 chemistry after the standard full-length cDNA (1–3 kb) library preparation. Three SMRT cells were sequenced on a PacBio Sequel system (Pacific Biosciences, USA). A SMRT Link version 6.0.0 was used for quality checking, filtering, and clustering of the Iso-Seq sequencing data.

A Flye assembler v2.9.1 (Kolmogorov et al., 2019) was utilized with the default parameters to create the de novo genome assembly of M. orangiensis. The Illumina data were used for polishing errors and improving the assembled genome by a Pilon v1.24 with the default parameters (Walker et al., 2014). HiRise software was used to conduct the chromosome-level genome assembly of M. orangiensis (Putnam et al., 2016). The assembled genome and the raw reads from the Hi-C libraries were aligned to construct a probability model for the genomic distance between read pairs. The genomic linking information between the contigs was generated, and mis-joins were corrected to create a pseudo-molecule-level scaffold genome. The Hi-C file comprising the contact matrices with duplication removal linking information was created using Juicer v1.6 (Durand et al., 2016b), and the Hi-C raw sequence data were aligned using BWA-MEM (Li, 2013). The contact map plot was drawn and visualized using Juicebox v1.6 (Durand et al., 2016a). Although there was no karyotype data from the Muraenolepididae family, we estimated that the scaffolds larger than 10 Mb were pseudochromosomes. The genome’s completeness was confirmed using Benchmarking Universal Single-Copy Orthologs (BUSCO) v5.4.6 (Manni et al., 2021) with the Actinopterygii odb10 database. The statistics for contig-level genome assembly were estimated by assemblathon_stats.pl (https://github.com/ucdavisbioinformatics/assemblathon2-analysis).

RepeatModeler2 (Flynn et al., 2020) was used to build the de novo repeat library with the default parameters in RepeatScout and RECON software. The consensus sequences and categorization information for each repeat were predicted using the Tandem Repeats Finder (Benson, 1999). All the repeats gathered by RepeatModeler were confirmed using the UniProt/SwissProt database (Dimmer et al., 2012). The raw LTR data from LTRharvest (Ellinghaus et al., 2008) and LTR_FINDER (Xu and Wang, 2007) were combined to construct the LTR library using LTR_retriever (Ou and Jiang, 2017). The highly accurate identification of long terminal repeat retrotransposons (LTR-RTs) was made using the LTR library. The repeating elements were identified using a de novo repeat library and RepeatMasker v4.1.4. The repetition landscape was estimated for each alignment using the Kimura distance.

The gene prediction of M. orangiensis was conducted using EvidenceModeler (EVM) v2.0 (Haas et al., 2008). First, the repeat masked genome was used for ab initio gene prediction using Augustus v3.5.0 (Stanke et al., 2006) and GeneMark-ES v4.68 (Lomsadze et al., 2005). Then, the hints for protein and ab initio prediction were extracted along with extensive protein sequences from the Actinopterygii clade of bony fishes in the UniProt/SwissProt protein database (UniProt, 2019) using ProtHint v2.6.0 (Bruna et al., 2020). The ProtHints output file was used to conduct the gene prediction for protein. Next, using the PASA pipeline v2.5.2 (Haas et al., 2008) with the Iso-Seq data, evidence for the transcriptome was generated for more accurate gene predictions of M. orangiensis. The results from the ab initio, transcriptome, and protein with their weights of contribution were integrated by EVM for the final gene predictions. Finally, to update the gene predictions into changing exons and adding un-translated regions (UTRs), we reused the PASA pipeline with the Iso-Seq data. Genome Annotation Generator v2.0.1 (Geib et al., 2018) was used to add the start and stop codon information and generate the gff3 file.

The predicted genes were annotated using NCBI BLAST v2.9.0 (Altschul et al., 1990) with a maximum e-value of 1e-5 by aligning the sequences to the NCBI non-redundant protein (nr) databases (Marchler-Bauer et al., 2011). InterProScan v92.0 (Paysan-Lafosse et al., 2022) was used with the default parameters and the protein sequences to identify the protein domains. Trinotate v3.2.2 (Bryant et al., 2017) was used for broad annotation of the transcriptome sequences, TransDeocder v5.5.0 was used to decode the amino acid sequences, SignalP v6.0 (Teufel et al., 2022) was used for the protein signal peptides, and Pfam (Sonnhammer et al., 1997) was used for the protein family predictions. The protein signatures were aligned against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to fetch KEGG-relevant functional annotations (Kanehisa and Goto, 2000). Gene Ontology (GO) terms were assigned to the genes for the gene function predictions and classifications using Blast2GO v6.0 (Conesa et al., 2005).

Nucmer in the MUMmer software package v4.02b (Marcais et al., 2018) with the following parameter -c 1100 -g 10000 and add –mum for unique match avoid repeat region was used to compare the genome sequences at the chromosome level. The long sequences corresponding to pseudochromosomes were selected and compared, excluding the unarranged scaffolds and contigs sequences. Circos (Krzywinski et al., 2009) was used for the comparative genomics, and a custom script was employed to convert the coordinate information from Nucmer into a Circos-compatible format. G. morhua was used as a reference genome, and the results of the chromosome comparison between the two genomes were diagrammed using Circos.

M. orangiensis and 20 well-assembled and annotated teleost genomes were utilized to identify the gene families and build a phylogenetic tree. The longest coding sequences (CDS) for each species were extracted from the annotated gff3 file using AGAT v1.0 (Dainat, 2021). To fetch the longest isoform, the agat_sp_keep_longest_isoform.pl script was used, and to extract the data to a fasta file, the agat_sp_extract_sequences.pl script was used. The CDS were converted to amino acid sequences using the “transeq” function in EMBOSS v6.6.0 (Rice et al., 2000). The protein sequences were used to classify the orthologous gene cluster using the OrthoMCL pipeline (Li et al., 2003) and applied the Markov Clustering Algorithm with default options in all the steps.

A phylogenetic tree was built using single-copy orthologous gene clusters. They were aligned using a Probabilistic Alignment Kit (Löytynoja and Goldman, 2005) with the codon alignment option. For the ensuing steps Gblocks (Castresana, 2000) were utilized to filter out the erratically aligned and gapped regions. Molecular Evolutionary Genetics Analysis (MEGA) v11.0 (Tamura et al., 2021) was used to construct the phylogenetic tree of 21 teleost species. The maximum-likelihood tree method was applied using RAxML (Stamatakis, 2014) under 1,000 bootstrap replications of the Jones-Taylor-Thornton (JTT) model. The divergence times calibration was conducted using Timetree (Hedges et al., 2006), and four teleosts were used (the minimum and maximum divergence times for Acanthochromis polyacanthus and Amphiprion percula were 27.0 and 29.3 million years ago (Mya); the smallest and largest divergence periods for G. morhua and Danio rerio were 180.0 and 264.0 Mya.

The previously identified orthologous gene clusters were utilized in conjunction with the five well-assembled and annotated Antarctic fish species to identify unique orthologs that live in frigid environments. Among the orthologous gene clusters, the specific gene families included in M. orangiensis were selected, and the functional categories were identified using a Blast2GO enrichment test with a Fisher’s exact test cutoff of a p-value of less than or equal to 0.05. CAFÉ v5.0 (Mendes et al., 2020) was used to perform the statistical tests for gene family expansion and contraction to identify how the gene families have altered during the evolution of M. orangiensis. The default parameters were used and identified the likelihood of the gene family expansion and contraction with a p-value of less than or equal to 0.05. The expanded and contracted genes of M. orangiensis were sorted, and the GO categories were identified using a Blast2GO enrichment test with a Fisher’s exact test cutoff of a p-value of less than or equal to 0.05.

M. orangiensis was classified as belonging to the Gadiformes order after comparing and analyzing the phylogenetic tree through the complete mitochondrial genome (Choi et al., 2020). Among the currently known teleosts in the Gadiformes order, it has been identified that there is a deficiency in MHC class II (Malmstrom et al., 2016). First, the presence of key genes in the cross-presentation pathway for the recognition of MHC class I, MHC class II, and pathogen-derived antigens was confirmed in the genome. Next, based on the analysis of the phylogenetic tree, the gene syntenies and chromosome locations of CD4 and CD74 genes, which are important MHC class II genes, were identified in Takifugu rubripes, G. morhua, and M. orangiensis. Although MHC class II is deficient in G. morhua, it has been specifically confirmed that the TLRs are repeatedly expanded (Sundaram et al., 2012). Based on the reference, the type and repetition of the TLRs in the genome were classified, and the location of the pseudochromosomes was identified. The TLRs in T. rubripes, G. morhua, D. mawsoni, and M. orangiensis were used for constructing the phylogenetic tree using MEGA-X based on the maximum-likelihood tree method with 1,000 bootstrap replicates. The TLR 8 gene from Homo sapiens was selected as an outgroup.

The read number and data size for the PacBio Sequel data were 398.66 Mb and 94.89 Gb. The Hi-C library had 676.79 Mb and 102.19 Gb for the number of reads and sequences. The number of reads and sequences of the Iso-Seq data were 0.18 and 38.66 Mb, respectively ( Supplementary Table 1 ).

This initial genome assembly was 890.95 Mb in length with a contig N50 of 0.73 Mb and the longest scaffold of 4.72 Mb, and the final chromosomal-level M. orangiensis genome assembly, which represented the first reference genome of family Muraenolepididae, has 893.75 Mb in length with a N50 scaffold length of 30.07 Mb and the longest scaffold of 39.77 Mb ( Table 1 ). Twenty-five high-quality pseudochromosomes were assembled and were listed in order from the longest to the shortest length. The longest pseudochromosome length was 39.77 Mb, and the shortest pseudochromosome length was 11.90 Mb ( Supplementary Table 2 ). The total length of the 25 pseudochromosomes was 83.5% of the total length of all genomic sequences, and the unassigned length was 147.62 Mb. The Hi-C interaction heat map was drawn and visualized using Juicebox ( Figure 1A ). BUSCO core genes were found as full BUSCO profiles in 3,398 (93.4%) of the 3,640 total BUSCO groups searched. Of these, 3,316 (91.1%) were single-copy, and 82 (2.3%) were duplicated BUSCOs. The fragmented and missing BUSCOs were 49 (1.3%) and 193 (5.3%), respectively ( Table 2 ).

The repetitive element sequences (REs) in the M. orangiensis genome were identified as tandem repeats and transposable elements (TEs). M. orangiensis have a variety of transposable elements subfamilies in the genome, and in total, the genome contained 51.28% transposable elements ( Figure 1B , Supplementary Table 3 ). Among the REs, 14.93% were long terminal repeats (LTRs), 9.98% were DNA transposons, 8.50% were long interspersed elements (LINEs), and 0.72% were short interspersed elements (SINEs). Among the LTR elements, “gypsy” was the highest after “unknown” at 4.92%, and among the DNA transposons, “CMC-EnSpm” was the highest at 3.30%. Among the LINEs, “Rex-Babar” was the highest at 4.03%, while in the SINEs, “MIR” was the highest at 0.30%. The total interspersed repeats had 1,613,229 counts 39.7% were masked in the M. orangiensis genome, and 10.04% were masked for simple repeats ( Supplementary Table 3 ). The Kimura distance was used to examine the TEs divergence distribution. The rates of the transitions and transversions were estimated on the alignment and converted to the Kimura distance to assess the age and transposition history of the Tes (Chalopin et al., 2015) with Kimura divergence values of 3 or below ( Figure 1B ).

A total of 34,553 genes in M. orangiensis were annotated, and the total length was 291.39 Mb. A total of 228,039 CDS and 229,440 exons were structurally annotated, and the lengths were 42.15 and 45.27 Mb, respectively ( Supplementary Table 4 ). Consequently, 34,854 protein genes were functionally annotated using the NCBI BLAST. Of the 34,854 protein genes, 27,691 (79.45%) were annotated in the GO database, whereas 7,345 (21.07%) were annotated in the KEGG database ( Supplementary Table 4 , Supplementary Figure 1A ). The Uniprot/Swiss-prot findings were recognized as 10,392 (29.81%), whereas Pfam and SignalP results were identified as 10,434 (29.94%) and 28,406 (81.5%), respectively ( Supplementary Table 4 ). The distribution of the top 20 GO terms in level 3 was identified. Among the biological process (BP) GO terms, organic substance metabolic process (GO:0071704) was the highest, followed by primary metabolic process (GO:0044238), and in molecular functions (MF), transferase activity (GO:0016740) was the highest, followed by hydrolase activity (GO:0016787). Organelle (GO:0043226) was the most common cellular component (CC) GO term, followed by intracellular anatomical structure (GO:0005622; Supplementary Table 5 , Figure 1C ). The enzyme codes for the genes were evaluated based on the GO annotation results, and 2,733 genes were matched with the transferases enzyme, which was the highest, followed by 2,461 genes matched with the hydrolases enzyme ( Supplementary Figure 1B ). The top species for the blast hits were G. morhua and Myripristis murdjan ( Supplementary Figure 2 ).

Utilizing 34,853 protein genes, the top 10 domains, families, sites, and repeats were analyzed using InterProScan. Among the top 10 domains, reverse transcriptase domain (IPR000477) was the highest, followed by NACHT nucleoside triphosphatase (IPR007111). The P-loop containing nucleoside triphosphate hydrolase (IPR027417) has the largest distribution in the top 10 family distribution, followed by immunoglobulin-like fold (IPR013783). In the top 10 repeats and site distribution, WD40 repeat (IPR001680) and protein kinase, ATP binding site (IPR017441) were the highest ( Supplementary Figure 3 ).

No known chromosomal-level assembly genomes exist for any Muraenolepididae family species; the homology between chromosomes was identified and studied using G. morhua in the adjacent Gadiformes order. The 23 chromosomes of G. morhua and the 25 pseudochromosomes of M. orangiensis were compared for the similarity of each chromosome to confirm the chromosomal stability of the Hi-C assembly. It was found that the portion of M. orangiensis corresponding to scaffold 197 did not exactly match the chromosomes of G. morhua and partially matched other chromosomes. The scaffolds 196 and 76 were matched in the form of splitting and dispersing on the G. morhua chromosomes ( Figure 2A ). Except for scaffolds 76, 196, and 197, most of the scaffolds showed high homology between the two species.

There were 4,519 multiple-copy orthologs, 2,653 unique orthologs, 17,607 other orthologs, and 8,598 unclustered genes in the M. orangiensis genome. The unique orthologs and unclustered genes outperformed the other fish species ( Figure 2B , Supplementary Table 6 ). The phylogenetic tree was built using 1,173 single-copy orthologs. The divergence time calibration was estimated using four fish species based on the phylogenetic tree, and the M. orangiensis divergence time was calculated to be 84.9 Mya. While no other Muraenolepididae species’ whole genomes have been sequenced, the phylogenetic tree revealed that M. orangiensis was connected to G. morhua in the Gadiformes order ( Figure 2C ).

M. orangiensis and four other Antarctic fish were compared for orthologous gene clusters ( Supplementary Table 6 ). There were 10,846 orthologous gene clusters identified to be common among the five Antarctic fish species, with M. orangiensis species alone detecting 624 orthologous gene clusters ( Figure 2D ). An enrichment test was conducted to identify the GO terms for the unique gene clusters identified in M. orangiensis, and as a result, five, including defense response to other organisms (GO:0098542), were identified as the highest in BP GO categories. Binding (GO:0005488) had the highest number of gene counts, with 353 in the MF GO category, followed by ion binding (GO:0043167; Supplementary Table 7 ).

Along the M. orangiensis lineage, 489 gene families that were significantly expanded and 741 gene families that were significantly contracted were identified. In comparison with the other Antarctic fish species included in the phylogenetic tree, M. orangiensis had less expanded gene families and more contracted gene families. At levels 5 and 6, the GO terms of the expanded and contracted gene families were analyzed using a gene enrichment analysis.

The highest BP in level 5 that related to the expanded gene families was organelle organization (GO:0006996) with 79 genes, followed by organonitrogen compound biosynthetic process (GO:1901566) and cellular macromolecule biosynthetic process (GO:0034645) both with 40 genes. The highest MF was DNA binding (GO:0003677) with 138 genes. For the CC, intracellular membrane-bounded organelle (GO:0043231) was the highest, with 355 genes. The highest BP and CC in level 6 were cytoskeleton organization (GO:0007010) with 66 genes and nucleus (GO:0005634) with 222 genes, respectively ( Supplementary Table 8 ). Among the contracted gene families, the highest BPs in level 5 were macromolecule modification (GO:0043412) and protein modification process (GO:0036211) with 32 genes, for the MF and CC, DNA binding (GO:0003677) and intracellular membrane-bounded organelle (GO:0043231) were the highest, respectively. The highest BP, MF, and CC in level 6 were RNA metabolic process (GO:0016070), pyrophosphatase activity (GO:0016462), and nucleus (GO:0005634), respectively ( Supplementary Table 9 ). The GO terms for the gene families that were quickly expanding and contracting covered a broad range of biological phenomena.

Although the functionality of the adaptive immune system was assumed to remain preserved since the origin of jawed vertebrates, substantial changes in the immune gene repertoire have lately been identified in teleosts (Malmstrom et al., 2016). The Atlantic cod has undergone one of the most drastic alterations, resulting in the full loss of the MHC class II pathway, which is important for detecting bacterial infections in vertebrates. Furthermore, this loss is accompanied by a significantly expanded repertoire of MHC I genes, which ordinarily produce molecules for defense against viral pathogens (Malmstrom et al., 2016). M. orangiensis, which likewise belongs to the Gadiformes order, was studied in terms of MHC class II-related immune systems.

The immune gene repertoires of the M. orangiensis genome were identified based on the findings of the annotation data and BLAST searches to identify the cause and scope of the MHC class II pathway gene loss seen in the Gadiformes order. The presence of 26 genes known to play important roles in MHC class I, MHC class II, and cross-presentation pathways was specifically looked for and compared with three other teleosts (Neefjes et al., 2011; Roche and Furuta, 2015). Initially, all the critical genes for MHC class I were found, with 23 genes identified, including 10 MHC I genes, which were the most identified in the M. orangiensis genome. In the three other teleost species, the MHC class I genes were analyzed, and it was confirmed that 73 MHC I genes were significantly duplicated in G. morhua. Seven antigen peptide transporter 2 (TAP2) genes were the second most identified MHC class I genes in M. orangiensis. Cross-presentation refers to the processing of foreign antigens into the MHC class I pathway (Heath and Carbone, 2001). A total of 11 important genes in cross-presentation were found, with two proteins identified: protein transport protein Sec61 subunit alpha (Sec61A), protein transport protein Sec61 subunit gamma (Sec61G), and cathepsin S precursor (CTSS). In addition, the signal sequence receptor subunit gamma (SSRG) gene is the sole gene missing from the cross-presentation key genes set. Lastly, except for two genes for endoplasmic reticulum aminopeptidase 1 (ERAP1) and one gene for autoimmune regulator (AIRE), no major MHC class II genes have been found in M. orangiensis. In contrast, in D. mawsoni and T. rubripes, CD4, CD74, and MHC class II, which are important key genes for MHC class II immunity, have all been identified. Given the loss of the CD4 and CD74 genes in M. orangiensis and G. morhua, which are critical for MHC class II immunization, these findings challenge the fundamental ideas regarding the evolution of the adaptive immune system and its components in vertebrates ( Table 3 ).

According to a previous study, the CD4 gene on chromosome 11 was fragmented, and the CD74 gene was lacking in the genome of G. morhua. We looked for a segment in M. orangiensis using the T. rubripes reference genome, which is known to have well-preserved CD4 and CD74 genes (Star et al., 2011). We identified well-preserved genes at both ends of the CD4 region, but the CD4 gene had been deleted. A conserved gene synteny of gene order was observed on pseudochromosome 10 surrounding the CD4 genes ( Figure 3A , Supplementary Table 10 ). The same method was used to identify the CD74 gene, which was not identified in the M. orangiensis genome, even among the well-preserved genes at both ends of the CD74 gene on pseudochromosome 10 ( Figure 3B , Supplementary Table 10 ). As a result, CD4 and CD74, the key genes of MHC class II, were not identified anywhere in the M. orangiensis genome.

The results of the CD4 and CD74 gene deletions in G. morhua have been reported previously, and one peculiarity was that TLRs were repeatedly expanded in the genome (Star et al., 2011). Accordingly, we identified TLRs in the M. orangiensis genome and compared them with three other teleosts. Notably, the TLR13 genes were highly expanded. A total of 44 TLRs were identified, with TLR13 having the highest gene counts, which were evenly distributed across six distinct pseudochromosomes, followed by seven TLR8 genes, which were identified on pseudochromosome 17 with the TLR7 gene ( Supplementary Table 11 ). Thirty-seven TLR13 genes were identified in G. morhua; however, only seven and five TLR13 genes were identified in D. mawsoni and T. rubripes, respectively. The TLR7 and TLR8 genes were not identified in the D. mawsoni and T. rubripes genomes. Finally, two TLR2 genes and one TLR3 gene were identified on pseudochromosomes 5, 10, and 24 ( Figure 4A ). Five different types of TLRs were used to construct a phylogenetic tree across three different fish species (G. morhua, D. mawsoni, and T. rubripes), with the TLR8 (AAF78036.1) gene from H. sapiens serving as an outgroup. Most of the TLRs were well-conserved, while TLR13 genes were significantly expanded in M. orangiensis, as they were in G. morhua ( Figure 4B ). There is no detailed study of this result in M. orangiensis, but compared with G. morhua, the absence of MHC class II, the immune system is assumed to have developed comprehensive mechanisms in both adaptive and innate immunity (Star et al., 2011).