An adult female specimen of M. coruscus (133.73 g) was collected in November 2018 from the Shengsi Islands in Zhejiang province, China. Gills were dissected and stored in liquid nitrogen until DNA extraction. Genomic DNA was extracted as previously reported (Venier and Canova, 1996) with small modifications. Briefly, gill samples were digested with RNase and proteinase K, adjusted to 2% SDS, and heated at 60°C for 10 min. Sodium perchlorate was added, and DNA was extracted once by gently shaking with 24:1 chloroform:isoamyl alcohol for 30 min at room temperature. DNA was then precipitated with cold 75% ethanol and suspended in TE buffer (10 mM Tris–HCl, 1 mM EDTA, pH 7.4). The quality and concentration of the extracted genomic DNA were checked using 0.3% agarose gel electrophoresis and a Qubit fluorimeter (Invitrogen, Carlsbad, CA, USA).

High-quality DNA was used for subsequent library preparation and sequencing using PromethION and Illumina platforms (Biomarker Technologies Corporation, Beijing, China). To obtain long non-fragmented sequence reads, ~15 μg of genomic DNA was sheared and size-selected (30–80 kb) with a BluePippin (Sage Science, Beverly, MA, USA). The selected fragments were processed using the Ligation Sequencing 1D Kit (Oxford Nanopore, Oxford, UK) according to the manufacturer's instructions and sequenced using the PromethION DNA sequencer (Oxford Nanopore, Oxford, UK) for 48 h. For the estimation and correction of genome assembly, an Illumina DNA paired-end (PE) library with an insert size of 350 bp was constructed in accordance with the manufacturer's protocol and sequenced with an Illumina HiSeq X Ten platform (Illumina, Inc., San Diego, CA, USA) with paired-end 150 read layout.

Reads from the two types of sequencing platforms were used at different assembly stages (Figure 1A). Long reads were filtered for length (>15,000 nt) and complexity (entropy over 15), while all short reads were filtered for quality (QC > 25), length (150 nt), absence of primers/adaptors, and complexity (entropy over 15) by using fastp (Chen et al., 2018). Using Jellyfish (Marçais and Kingsford, 2011), the frequency of 31-mers in the Illumina filtered data was calculated with a 1 bp sliding window (Vurture et al., 2017). Long reads were then assembled using wtdbg2 (Ruan and Li, 2020), which uses a fuzzy Bruijn graph. As it assembles raw reads without error correction and then builds the consensus from intermediate assembly output, multiple stages of error correction, gap closing, and polishing were applied. Original output was realigned against the long reads and polished using Minimap2 (Li, 2018) and Racon (Vaser et al., 2017), first with filtered reads to bridge potential gaps, and then with the filtered reads to correct for error. Finally, Pilon (Walker et al., 2014) was used to polish and correct for sequencing error using the short reads.

We used Braker (Hoff et al., 2019) to perform ab initio gene prediction, combining methods that integrate ab initio gene prediction and RNA-seq-based prediction to annotate the protein-coding genes in the M. coruscus genome. These raw RNA-seq reads were downloaded from the EBI for three independent transcriptomic projects, covering multiple tissues and multiple conditions: PRJNA301064 (Xu et al., 2016), PRJNA269003 & PRJNA269004. The resulting predictions were then filtered for the presence of at least one InterPro (Jones et al., 2014) pattern using InterProScan (Mitchell et al., 2019).

The transposable elements were annotated using a de novo prediction by using RepeatModeler (Smit and Hubley, 2017) and LTR-Finder (Stanke et al., 2008). The repetitive sequences returned from these two algorithms were combined to compile a non-redundant repeat sequence library. With this library, we scanned the representative sequences in the M. coruscus genome using RepeatMasker (Smit et al., 2019).

The completeness of gene regions was further assessed using BUSCO (Simão et al., 2015), using a Metazoa (release 10) benchmark of 954 conserved Metazoa genes.

To assess the macro-synteny between M. coruscus and M. galloprovincialis (Murgarella et al., 2016) genomes, we reciprocally mapped all M. coruscus and M. galloprovincialis scaffolds with a minimum overlap of 100 nt.

The versions, settings, and parameters of the software used in this work are as follows:

Genome assembly: (1) fastp: version 0.20.0, short-read parameters: -q 25 -y -Y 15 -l 150 -detect_adapter_for_pe; (2) fastp: version 0.20.0, long-read parameters: -Q -y -Y 15 -l 15000; (3) wtdbg2: version 2.4, parameters: -x rs -k 23 -p 0 -AS 6 -R -g 1567m -rescue-low-cov-edges; (4) wtpoa-cns: version 2.4, default parameters; (5) minimap2: version 2.17, parameters: -x map-ont -r2k; (6) racon: version 1.4.3, default parameters; (7) bwa: version 0.7.17, mode mem, default parameters; (8) pilon: version 1.23, parameters: -diploid -fix all -changes; (9) BUSCO: version 4.0.2, parameters: -l metazoa_odb10; (10) RepeatModeler: version 1.0.11, parameters: -database mussel; (11) LTR_Finder: version 1.07, default parameters; (12) RepteatMasker: version 4.0.9, parameters: -lib mussel-families.fa; (13) Braker: version 2.1.4, parameters: -gff3 -softmasking; (14) InterProScan: version 5.42-78.0, parameters: -f tsv -dp -iprlookup -goterms.

K-mer analysis: (1) jellyfish: version 2.3.0, parameters: -m 31 -C -s 10G; (2) GenomeScope: version 2.0, default parameters.

Mitochondria annotation: (1) MITOS: revision 999, online version, parameters: “Genetic code 5”.

Synteny blocks: (1) nucmer: version 3.23, parameters: -simplify -maxgap = 500 -mincluster = 100; (2) circos: version 0.69-9, default parameters (tutorial 5.9).

After sequencing with the PromethION platform, a total of 11.31 million (161.04 Gb) long reads were generated, and these were used for the subsequent genome assembly. The N50 of the sequences produced was 21,771 nt. The Illumina HiSeq X Ten platform produces 288.22 million (86.47 Gb) paired-ended short reads (150 nt). The genome size of M. coruscus was estimated to have 2n = 28 chromosomes (Pérez-García et al., 2014) and a C-value of 1.90 pg (Ieyama et al., 1994) or 1.85 Gb; therefore, the average sequencing coverage was 87× and 46×, respectively (Table 1).

The frequency of 31-mers in the Illumina filtered data was calculated using Jellyfish and followed a Poisson distribution (Figure 1B). The proportion of heterozygosity in the M. coruscus genome was evaluated to be 1.64%, and the genome size was estimated as 1.57 Gb, with a repeat content of 36.35% (Table 2). Long-read assembly using wtdbg2 (Ruan and Li, 2020), polished using Racon and sequence-corrected using short reads and Pilon, produced an assembled genome of M. coruscus containing 10,484 contigs with a total length and contig N50 of 1.90 Gb and 898 kb, respectively (Table 2).

The transposable elements and repetitive sequences were annotated using RepeatMasker and LTR-Finder. In total, 1.01 Gb (52.83%) of the genome was identified as repetitive sequences (Table 3).

We used a combined method that integrates ab initio gene prediction and RNA-seq-based prediction to annotate the protein-coding genes in the M. coruscus genome. In total, 42,684 distinct gene models were annotated.

To estimate the quality of the hybrid genome assembly, short reads were aligned to the consensus genome, and 98.42% did align overall, suggesting that our assembly results contained comprehensive genomic information.

The completeness of the gene models was also assessed using BUSCO (Simão et al., 2015) using a Metazoa (release 10) benchmark of 954 conserved Metazoa genes: 91.09% had complete gene coverage (including 13.63% duplicated ones), 5.35% were fragmented, and only 3.56% were missing (Figure 1C). This largely supports a high-quality M. coruscus genome assembly and gene models.

To uncover the macro-synteny between the M. coruscus and M. galloprovincialis (Murgarella et al., 2016) genomes, we reciprocally mapped all M. coruscus and M. galloprovincialis scaffolds with a minimum overlap of 100 nt (Table 4). A total of 28.12% of M. galloprovincialis scaffolds mapped to a M. coruscus scaffold while 70.25% M. coruscus scaffolds mapped to a M. galloprovincialis scaffold. Figure 2A reports the synteny blocks of over 3,000 nt; the map shows a good synteny between the two genomes despite M. galloprovincialis genome fragmentation.

The mitochondrial genome was manually recovered from the genome assembly. The contig was validated for continuity and circularity and annotated using MITOS (Bernt et al., 2013). The complete mitochondrial genome (Figure 2B) was compared to the reference M. coruscus genome (Lee and Lee, 2016). Only one haplotype was recovered, which differed from the reference by only 6 SNPs and one 2-nucleotide insertion.