The tissue sample was obtained from a single wild E. affinis female caught off the coast of Ainan, Ehime Prefecture, Japan. High-molecular-weight genomic DNA was isolated from fin using NucleoBond^® AXG columns with NucleoBond^® Buffer Set IV (Macherey-Nagel, Düren, Germany). The quantification of gDNA was performed by Quant-iT™ dsDNA Broad-Range Assay Kit (Invitrogen, Carlsbad, CA, USA), and molecular weight was estimated on 0.75% agarose gel by pulsed-field electrophoresis. The whole-genome sequencing library was prepared using Chromium Genome Library & Gel Bead Kit v. 2 (10x Genomics, Pleasanton, CA, USA) as described in the Chromium Genome Reagent Kit v. 2 User Guide. The library was sequenced on Illumina HiSeq X sequencing system using pair-end (2 × 150 bp) sequencing. The sequencing generated 920.3 million (M) reads of total 138.04 Gb with 93.8% and 87.7% of base having quality score Q > 30 in R1 and R2 reads, respectively.

Oxford Nanopore Technology (ONT) sequencing was used to obtain long reads for scaffolding. Libraries were generated using standard protocols from ONT with the SQK-LSK108 ligation sequencing kit (Oxford Nanopore Technologies, Oxford, United Kingdom). GridION X5 sequencing was performed according to the manufacturer's guidelines using three independent FLO-MIN107 (R9.5) flow cells. Base-calling was done by Albacore v. 1.2.4. (Oxford Nanopore Technologies, Oxford, United Kingdom). Raw reads were filtered by quality value QV9 (--minqual 9), and heads and tails were 50 bp trimmed each (--headtrim 50; --tailtrim50) in Yanagiba v. 1.0.0 (Taranto, 2017). The ONT sequencing generated 1.56 M reads with a total of 15.02 Gb, an average length of 9,650 bp, and N50 = 18,747 bp.

Prior to genome assembly, genome characteristics were estimated based on Illumina reads. Jellyfish v. 2.2.6 (Marçais and Kingsford, 2011) was applied to generate k-mer counting and frequency distributions of 19-, 21-, and 23-mers. Genome size, heterozygosity, and repeat content were estimated based on the generated k-mer count distributions using GenomeScope (Vurture et al., 2017) with high frequency k-mer cutoff = 10,000. The estimate in GenomeScope is based on an equation that models four evenly spaced negative binomial distributions of the k-mer profile to measure the relative abundances of heterozygous and homozygous, unique, and two-copy sequences (Vurture et al., 2017). The genome size estimate ranged from 745.82 Mb (k = 19) to 755.80 Mb (k = 23), heterozygosity rate was roughly estimated to be 0.67% (67 SNPs per 10 Kb), and repeat content estimate ranged from 137.00 Mb (k = 23) to 170.35 (k = 19; Supplementary Table 1).

Assembly of Illumina reads was performed in Supernova assembler v. 1.1.5 (10x Genomics, Pleasanton, CA, USA), with default parameters, except maximum reads (--maxreads), set at 386 M input reads to achieve 56 × raw coverage, as suggested in the Supernova protocol. Assembled scaffolds were loaded together with filtered ONT reads into PBJelly v. 15.8.24 (English et al., 2012) where identification of gaps >25 bp, gap filling, and scaffolding were performed using default parameters. In total, 11.31 Mb of gaps (34.8% of initial gap size) were successfully closed. The 787.60-Mb initial draft genome assembly was presented in pseudohaplotype format and consisted of 19,850 scaffold sequences (>1 Kb), of which 21.17 Mb (2.69%) represented unknown bases.

The redundancy of the genome assembly was reduced in three steps: First, 8,191 duplicated scaffolds (24.24 Mb) were removed using dedupe.sh from BBTools v. 38.87 (Bushnell, 2014). Furthermore, scaffolds < 2 Mb were clustered in CD-HIT v. 4.8.1 (Li and Godzik, 2006) with identity threshold ≥ 99% (-c 0.99) and word size -n 10. This step removed 1,044 scaffolds (5.1 Mb). Finally, retained scaffolds were self-aligned in LastZ v. 1.04 (Harris, 2007) with alignment identity threshold ≥ 99% (--identity = 99) and query coverage threshold ≥ 95% (--coverage = 95). If two different scaffolds were self-aligned, the longer one was retained and the shorter one was discarded. In total, 9,250 scaffolds (29.4 Mb; 3.7%) were removed from the initial assembly because of potential duplication or redundancy. The final scaffold-level assembly consisted of 10,600 scaf-folds of 758.20 Mb (19.89 Mb gaps), N50 = 24.78 Mb, and the longest scaffold = 35.13 Mb.

We obtained two types of linkage evidence: diploid, based on full-sib family linkage analyses, and haploid, based on linkage analyses of interspecific hybrids (Yoshitake et al., 2018) of E. affinis female and T. orientalis male. Based on a genome coordinate of each marker in the linkage maps, we anchored and oriented scaffolds into pseudochromosomes.

To construct a diploid linkage map, DNA was extracted from fin clips of parents and 94 of their progeny using NucleoSpin^® Tissue (Macherey-Nagel, Düren, Germany). All specimens were assessed for body weight, standard length, head length, and body depth (Supplementary Table 2). Genotyping by random amplicon sequencing-direct (GRAS-Di^®) libraries of each specimen were prepared according to the protocol of Hosoya et al. (2019). The final PCR products were pooled, purified using the MiniElute PCR Purification Kit (Qiagen, Hilden, Germany), and applied for pair-end (2 × 76 bp) sequencing on the Illumina NextSeq 500. Library preparation and sequencing were done by Bioengineering Lab. Co., Ltd., under a license agreement, as GRAS-Di^® is patented by the Toyota Motor Corporation (Aichi, Japan) (patent ID P2018-42548A). The sequencing generated ~160 M (12.17 Gb) raw reads. These were trimmed in TrimGalore v. 0.6.4 (Krueger, 2019) to remove residues of indexes, nucleotides with quality ≤ Q20, and reads ≤ 20 bp after trimming. In total, 64.81 M reads in pairs (675,154 reads per sample on average) with a total length of 9.67 Gb (100.72 Mb per sample on average) were retained after trimming (Supplementary Table 2). Trimmed reads of each sample were mapped onto reference assembly using Bowtie2 v. 2.4 (Langmead and Salzberg, 2012) allowing no mismatches in the read seed (-N 0). The overall mapping rate was 92.2% with 74.9% of reads mapped only once (Supplementary Table 2). Sorted BAM files were created using samtools v. 1.10 (Danecek et al., 2021). Variants were called in bcftools v. 1.10.2 (Danecek et al., 2021) by mpileup and call commands using multiallelic-caller and variants-only flags (-mv). Filtering of SNPs was performed in vcftools v. 0.1.16 (Danecek et al., 2011) removing indels (--remove-indels), SNPs with quality ≤ 30 (--minQ 30), SNPs not in HWE (--hwe 0.05), and genotypes with depth ≤ 10 (--minDP 10). Minor allele frequency was set to 0.05 (--maf 0.05), retaining 5,853 SNPs of 368,265 raw variants. Genotypes were phased, and both female and male linkage maps were constructed in TMap v. 1.1 (Cartwright et al., 2007). All segregating markers that showed polymorphism in at least one parent were used. The ratio of marker segregation was calculated by chi-squared test. Markers showing significantly distorted segregation (p-value < 0.001) were excluded from the map construction. A minimum logarithm of odds (LOD) threshold of 5.0 was selected to assign markers to 24 linkage groups. Recombination rates were calculated by the multipoint-likelihood maximization, and map distances were converted by Kosambi mapping function. A total of 852 of 1,332 polymorphic loci were assigned to 24 linkage groups covering a total length of 1,554.7 cM of the E. affinis genome (Supplementary File 1).

To obtain haploid linkage evidence, interspecific hybrids of E. affinis and T. orientalis were produced. Briefly, eggs of a single E. affinis female were mixed with cryopreserved sperm of T. orientalis provided by the Aquaculture Research Institute, Kindai University, Japan. Seawater was immediately added to activate gametes and induce fertilization. After 1–2 h, eggs with proceeding cleavage were selected and transferred to the hatching tank at 24°C. After 36–40 h, 202 individuals that hatched or died after reaching somite formation were separately transferred to a 1.5-ml tube and stored in 100% ethanol. DNA was extracted from parents and F1 hybrids using a NucleoSpin Tissue XS kit (Macherey-Nagel, Düren, Germany). A sequencing library was prepared from 135 specimens, which provided sufficient DNA for use of the Nextera DNA Library Preparation Kit and Nextera Index Kit (Illumina) following the manufacturer's protocols. The library was sequenced on two lines of Illumina HiSeq X sequencing system. The sequencing resulted in 2,274 M sequencing reads with a total of 342 Gb. Reads were mapped onto a reference obtained by combining E. affinis scaffold-level assembly and T. orientalis genome (Suda et al., 2019) using BWA mem v. 0.7.15 (Li and Durbin, 2009). Mapped reads were sorted in samtools v. 1.10 (Danecek et al., 2021), and variants were called in bcftools v. 1.10.2 (Danecek et al., 2021) by mpileup and call commands. Linkage evidence of scaffolds was obtained through linkage analysis of hybrids in SELDLA v. 2.0.9 (Yoshitake et al., 2018) using 13,403,357 SNPs specific to E. affinis.

Female and male linkage maps and scaffold linkage evidence were transformed to BED files and merged. Pseudochromosomes were then reconstructed using ALLMAPS v. 1.1.7 from the JCVI utility libraries v. 0.7.5 (Tang et al., 2015) with inter-scaffold gaps set to a fixed size of 100 Ns. The package was used to merge bad files and to anchor, order, and orient genomic scaffolds using default parameters. Overall, 387 scaffolds with total length of 685.79 Mb (90.7% of scaffold-level assembly) were anchored onto 24 pseudochromosomes leaving 10,213 scaffolds of total length 72.42 Mb unplaced. Only two unplaced scaffolds had length > 1 Mb (Supplementary Figure 1). The final assembly contained 10,237 scaffolds of 758.24 Mb (19.96 Mb gaps), N50 = 29.18 Mb, and longest scaffold = 35.73 Mb (Table 1).

Total RNA was extracted from eight tissues (brain, liver, kidney, ovary, testis, spleen, gill, muscle, and intestine) using TRIzol (Invitrogen, Carlsbad, CA, USA) and the NucleoSpin^® RNA Plus extraction kit (Macherey-Nagel, Düren, Germany) following the manufacturer's protocols for each tissue. RNA extracts were quantified using a NanoPhotometer N50 (Implen, München, Germany) and subsequently combined in equimolar quantities into a single pool for sequencing. RNA sequencing library was prepared by MGIEasy RNA Directional Library Prep Set (MGI Tech Co Ltd.). Pair-end (2 × 150 bp) sequencing was performed on a DNBSEQ-G400 sequencer (MGI Tech Co Ltd.). All procedures were conducted according to the manufacturer's protocols. The sequencing generated ~682 M (102.53 Gb) raw reads. Quality metrics for sequencing reads were initially examined in FastQC v. 0.11.9. (Andrews, 2020). Rare, possibly erroneous, k-mers were removed in Rcorrector v. 1.0.4 (Song and Florea, 2015) with default parameters, and adapters and low-quality bases were trimmed in TrimGalore v. 0.6.4 (Krueger, 2019) with parameters --length 36 -q 5 --stringency 3 -e 0.1 retaining ~614 M (92.18 Gb) pair-end reads. The FastQC results revealed deviation from normal distribution of GC content, and a high number of overrepresented sequences, possibly due to incomplete polyA capture during library preparation. Thus, trimmed reads were mapped against ribosomal RNA (rRNA) sequence database SILVA release 128 (Quast et al., 2012) using Bowtie2 v. 2.4 (Langmead and Salzberg, 2012) with parameters --nofw --quiet -D 20 -R 3 -N 0 -L 20 -i S,1,0.50 to remove rRNA contaminants. The 359.6 M (53.9 Gb) reads that did not map (--un-conc-gz) to the SILVA database were processed for de novo assembly in Trinity v. 2.11.0 (Grabherr et al., 2011) with default k-mer size 25 and --SS_lib_type RF --min_contig_length 300 flags. Initial transcriptome assembly resulted in 271,656 contigs of 242.43 Mb, N50 = 1.20 Kb. These were transferred to super transcripts (Davidson et al., 2017) by Trinity's v. 2.11.0 Trinity_gene_splice_modeler.py, and further redundancy was reduced by Bellerophon pipeline v. 1.0 (Kerkvliet et al., 2019), removing minimally expressed (transcripts per million cut off = 1) and highly identical (95%) contigs (CDHIT-EST -c 0.95). Final transcriptome assembly consisted of 49,510 contigs of 94.58 Mb, N50 = 3.62 Kb, with the longest contig = 56.88 Kb (Table 1).

A de novo repeat library was generated using RepeatModeler v. 2.0.1 (Flynn et al., 2020) and MITE Tracker v. 1.0.0 (Crescente et al., 2018) with default parameters. The genome was then screened for repeats and low complexity regions by RepeatMasker v. 4.1.1 (Smit et al., 2015) in two runs using (i) de novo-generated repeat library and (ii) a Dfam database of interspersed repeats, release 3.3 (Storer et al., 2021). Results of the runs were analyzed together to generate final non-redundant repeat annotation. Repetitive regions accounted for 25.59% (194.01 Mb) of genome assembly (Figure 1A; Table 1). These included 12.42% unclassified repeats, 3.95% retrotransposons, 5.81% DNA transposons, 0.08% small RNAs, 0.01% satellites, 2.61% simple repeats, and 0.42% low complexity regions (Supplementary File 2).

Gene models were predicted in MAKER v. 3.01.03 (Holt and Yandell, 2011) in three successive runs. Prior to the first run, complex repeats were retrieved from the repeat annotation file and submitted to MAKER as pre-identified repeat elements (rm_gff) while still enabling the software to identify and soft mask simple repeats internally (Card, 2017). In this approach, complex repeats are hard masked so that they do not confound the ability to identify coding genes, while simple repeats remain available for inclusion in gene annotations, as many protein-coding genes contain runs of low-complexity sequence (Toll-Riera et al., 2011). During the first run, the E. affinis transcripts were aligned to the genome by BLASTN (Camacho et al., 2009) and protein sequences of Danio rerio, Gasterosteus aculeatus, Hippocampus comes, Oreochromis niloticus, Oryzias latipes, Seriola dumerili, Sparus aurata, and Takifugu rubripes from the Ensembl database v. 103 (Yates et al., 2019) along with Thunnus orientalis (Yasuike et al., 2016) by BLASTX (Camacho et al., 2009). Subsequently, BLAST hits were polished by Exonerate v. 2.4.7. (Slater and Birney, 2005) est2genome and protein2genome. All filtering statistics for BLAST and Exonerate were as the default by MAKER. The second and third runs of MAKER utilized gene models from the first, followed by the second, runs to train ab initio gene prediction tools SNAP v. 2013-02-16 (Korf, 2004) and Augustus v. 3.3.3. (Stanke and Waack, 2003). This bootstrap process allows to iteratively improve the performance of ab initio gene predictors as they require existing gene models on which to base prediction parameters. SNAP was retrained using gene models with an annotation edit distance (Holt and Yandell, 2011) (AED) ≤ 0.25 and amino acid length of ≥ 50. BUSCO v. 5.1.2 (Simao et al., 2015) with --long argument and actinopterygii_obd10 lineage dataset was used to retrain Augustus using genomic regions of RNA annotations from the previous run including an additional 1,000 bp on each side as input file. Both SNAP and Augustus were run with default parameters specified in MAKER. Only gene models with AED <0.5 were retained in the final annotation set.

For functional annotation of predicted genes, predicted protein sequences were mapped against UniProtKB/Swiss-Prot (The UniProt Consortium, 2021) and NCBI non-redundant (O'Leary et al., 2016) protein databases using BLASTP (Camacho et al., 2009) with an e-value threshold of 1e−6. Additionally, protein motifs, domains, and signatures were annotated using Interproscan v. 5.48 (Jones et al., 2014), and Gene Ontology (GO) terms were obtained from the corresponding InterPro entry. Kyoto Encyclopedia of Genes and Genomes Orthologs (KOs) were assigned to predicted proteins using KofamKOALA (Aramaki et al., 2019) with an e-value threshold of 1e−3.

In total, 23,059 putative genes spanning 32.0% of the genome were predicted (Figure 1A; Table 1). We found 21,313 (92.4%) predicted genes to match at least one of the databases (Table 1), and at least one GO and/or KO term was retrieved for 11,429 and 14,796 predicted genes, respectively.