The family Scombridae is divided into 2 subfamilies (Gasterochismatinae and Scombrinae), with 15 genera and around 49 described species, comprising mackerels, bonitos, and tunas [1]. The representative genus of the Scombridae, i.e., Scomber, includes 4 species: S. scombrus, S. japonicus, S. australasicus, and S. colias. The Atlantic chub mackerel, Scomber colias (Gmelin, 1789) (NCBI:txid338315, FishBase ID:54736), is a small coastal pelagic fish that is distributed widely, being found in the Atlantic Ocean from the Bay of Biscay to South Africa (including the Canary, Madeira, Azores, and Saint Helena Islands), and in the Mediterranean Sea (Figure 1) [2]. Scomber colias is usually found at depths of up to 300 m and occupies a key position in the trophic web. This species acts as a link between primary producers and top predators, since it feeds mainly on zooplankton and small pelagic fish, and is an essential element of the diet of larger pelagic fish (e.g., tuna, swordfish, and sharks) and marine mammals (e.g., dolphins and seals) [3]. Besides its ecological importance, S. colias also supports important commercial fisheries for several countries across its distribution range, being an important component in the diet of local populations [1, 4]. This is probably related to its nutritional value, as this mackerel is a valued source of important fatty acids for human nutrition, particularly docosahexaenoic acid (DHA), an omega-3 fatty acid [5, 6]. Additionally, S. colias is used as bait for the tuna longline and handline fisheries, and is caught in purse seine and pelagic trawl fisheries which target sardines and anchovies [7].

The availability of S. colias makes it a sustainable marine resource [6] and a viable alternative to the European sardine (Sardina pilchardus), which is under fishing restrictions due to a population decline. Curiously, fluctuations in abundance and a northwards shift in the distribution of S. colias, with a likely inverse relationship with sardine abundance, have been recently demonstrated [8]. Due to its ecological and economic importance, S. colias has been the focus of several recent studies on different aspects of its fisheries and biology [3, 8, 9]. Yet, genomic resources for the species are still limited. Presently a (liver) transcriptome [10], a mitogenome  [11], and single-nucleotide polymorphism (SNP) data obtained through restriction site-associated DNA sequencing [12], have been described for the species. With the vast majority of the world’s fish stocks already in collapse, and with climate change as additional pressure, information on fish genomes is becoming a pressing tool to address conservation efforts [13, 14]. Here, we report the first high-quality draft genome of S. colias, assembled with Illumina and Pacific Biosciences (PacBio) Single Molecule High-Fidelity (HiFi) reads. This resource provides a critical platform to uncover the species’ adaptive physiological potential in a changing environment. Specifically, it will help understand the current observed populational northward shift, postulated to be part of a more general expansion of species from warmer areas [8]. Moreover, being one of the genomes with higher quality within the family Scombridae and the first within the Scomber genus, this information will help to improve the conservation, management, and sustainable exploitation of this valuable fish resource as well as that of its highly valued congeners.

Two specimens of S. colias were collected at 2 sampling time points. The first specimen was collected in 2017, during the “Programa Nacional de Amostragem Biológica” managed by the Instituto Português do Mar e da Atmosfera” (IPMA), in North Atlantic waters (41.501944 N 8.851667 W). From this individual, 2 tissue types were collected, and were stored in 100% ethanol (muscle) or RNA later (liver). Liver tissue was used to produce and describe the first liver transcriptome of S. colias [10]. Muscle tissue was used in the present study, for genomic DNA (gDNA) extraction using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. The gDNA was then used for Illumina paired-end (PE) sequencing (described below). The second specimen was caught in 2020, near Mira, Portugal (40.5588270 N 9.4529720 W). Immediately upon harvesting, the muscle was snap frozen in liquid nitrogen. The frozen tissue was shipped to Brigham Young University DNA Sequencing Center (BYU), where gDNA with high molecular weight was extracted from 1.1 g of muscle using the QIAGEN Genomic-tip 20/G kit. The quality and concentration of the gDNA were assessed with Qubit Fluorometric system (ThermoFisher), and the fragment size was determined with a fragment analyser (Agilent Technologies, RRID:SCR_013575) before loading on the Pacbio Sequel II system (PacBio Sequel II System, RRID:SCR_017990).

For the first DNA sample, Illumina PE library preparation and sequencing were carried out by Macrogen, Inc. (Seoul, Korea), using Illumina HiSeq X Ten platform (Illumina HiSeq X Ten, RRID:SCR_016385), with 250 bp PE configuration. For the second specimen, PacBio HiFi library preparation and sequencing were performed at BYU, following the manufacturer’s recommendations [15]. The size-selected fraction had a mean read length of 15.3 kbp and was selected on the SageELF system (Sage Science, RRID:SCR_014808). The sequencing was conducted on two single-molecule, real-time (SMRT) cells using Sequel II system v.9.0, with a run time of 30 h, and 2.9 h pre-extension. The circular consensus analysis was performed in SMRT® Link v9.0 [16] under default settings (the statistics of raw data generated from each PacBio SMRT cell can be viewed in Additional File 1 [17, 18]).

Both short- and long-read datasets were assessed by FastQC v.0.11.8 software (FastQC, RRID:SCR_014583). Trimmomatic v.0.38 software (Trimmomatic, RRID:SCR_011848) [19] was used to filter and remove low quality reads as well as the adaptors of the Illumina dataset (LEADING:5 TRAILING:5 SLIDINGWINDOW:4:20 MINLEN:50). Next, trimmed datasets were used to check the overall characteristics of the S. colias genome (i.e., genome size, heterozygosity, unique content), through GenomeScope 2.0 [20]. Briefly, Jellyfish v.2.2.10 software (Jellyfish, RRID:SCR_005491) [21] was used to build k-mer frequency distributions, and the final k-mer counts (k-mer 21, 25, 31) were submitted to the GenomeScope 2.0 online platform. On the other hand, HiFi reads were filtered in two ways (Figure 2). First, mitochondrial reads were removed by BLAST searches (BLASTN, RRID:SCR_001598) using a prebuilt database of mitochondrial sequences (database build protocol: (1) select all complete mitogenomes present in the nucleotide database of the National Center for Biotechnology Information (NT-NCBI); (2) select by taxon (Actinopterygii; txid:7898); (3) sequence length filter 15,000–50,000 bp; (4) build a database with the makeblastdb application of NCBI-BLAST+ v.2.9.0). Second, to filter out possible sources of contamination (artefactual or biological), HiFi reads were checked by BLAST (BLASTN) against NT-NCBI. Only HiFi reads with match hits over 90% identity and query coverage of 50% in the Actinopterygii taxon (NCBI:txid7898), or without match hits at all, were considered for further analysis (Figure 2).

Given that 2 specimens were used for the distinct sequencing approaches, i.e., PacBio HiFi and Illumina PE, the whole mitochondrial genome (mtDNA) was assembled and characterised for both specimens. For specimen 1, trimmed Illumina PE reads were used to assemble mtDNA in GetOrganelle v.1.7.1 [22] with optimised parameters (-F animal_mt -w 121 -R 10 -k 85,95,105,115,125) (Figure 2). For specimen 2, a new pipeline was designed to produce a mtDNA assembly from the PacBio HiFi long reads (Figure 2). The PacBio HiFi mtDNA reads, previously filtered (see above), were corrected using Hifiasm v.0.13-r308 (Hifiasm, RRID:SCR_021069) [23] with optimised parameters (–write-ec). Since Hifiasm is not optimised to assemble circular molecules (which are expected for mtDNA), the corrected PacBio HiFi mtDNA reads were assembled using Unicycler v.0.4.8 [24], a software package designed to assemble bacterial genomes and so optimised for circular assemblies, with default parameters. Annotation and visual representation of both mtDNA assemblies were produced using MitoZ v.2.3 [25] with optimised parameters (–genetic_code 2; –clade Chordata; –topology circular), using the PE reads for coverage plotting. Furthermore, annotations were manually validated by comparison with other mitochondrial genomes of the genus Scomber, available at NCBI (see Data Availability).

For whole-genome assembly a combined approach, using short- and long-read assemblies, was applied (Figure 2). While long-read assemblies were mainly used to produce the primary assembly, short-read assemblies were used to scaffold and improve the contiguity of the basal assembly. In summary, short-read assemblies were performed with the W2RAP pipeline v.0.1 [26], following the authors’ protocol. First the k-mer analyses toolkit (KAT) v.2.4.1 software (KAT, RRID:SCR_016741) [27] hist module was applied to determine the ideal k-mer cut-off, before W2RAP with optimised parameters (-t 30; -m 500; –min_freq 14; -d 32; –dump_all 1; -k: 144, 180, 200, 224) was used to produce 4 assemblies (Figure 2). To generate the long-read assembly, multiple software and parameters were initially tested. PacBio HiFi reads were assembled in Hifiasm v.0.13-r308 [23] with a range of parameters (k = 21, 25, 31, 41, 45, 51; l = 0, 2) and in HiCanu v.2.1.1 [28] with optimised parameters (default). While Hifiasm generated 2 pseudo-haplotypes per assembly, HiCanu generated 1 merged assembly. To choose the “best” assembly we applied a series of analyses, including Bandage (a bioinformatics application for nagivating de novo assembly graphs easily) v.0.8.1 [29] and manual inspection; Benchmarking Universal Single-Copy Orthologs (BUSCO) v.5.2.2 (BUSCO, RRID:SCR_015008) [30] with Eukaryota and Actinopterygii databases was used to assess the gene completeness of the assemblies, and Quality Assessment Tool for Genome Assemblies (QUAST) v.5.0.2 (QUAST, RRID:SCR_001228) [31], to check general metrics of the assemblies (Figure 2). Due to discrepancies in the length of the Hifiasm primary and alternative pseudo-haplotypes, we chose to concatenate them in a single assembly. At this point, the assembly with the highest complete BUSCO scores, highest contiguity (N50), and longest contig, was selected for further analysis. The pseudo-haplotypes were separated by purge_dups v.1.2.5 (purge dups, RRID:SCR_021173) [32]. After the first round of purging and inspection by k-mer plot, produced by the KAT tool, cutoffs were manually adjusted. To assess the influence of purge_dups in the genome, BUSCO (rate of deduplicates) and QUAST (N50 and genomic length per pseudo-haplotype) were used. Next, to improve the contiguity and quality of the assembly, short-read assemblies were used to structurally scaffold the assembly without the introduction of any new bases in the assembly, similar to the literature [33, 34] (Figure 2). The 4 short-read assemblies were inputted to the Long Interval Nucleotide K-mer Scaffolder (LINKS) v.1.8.7 [35], being used as long reads; using several distance values, i.e., -d 0.5, 1.5, 3, 9, 27 kb, the primary assembly was rescaffolded interactively for 5 rounds (additional parameters: -k 21 -e 0.5). Furthermore, the scaffolded genome and the long-read assemblies, initially produced by Hifiasm and HiCanu and discarded based on contiguity and completeness, were inputted to Cobbler v.0.6.1 [36] and RAILS v.1.5.1 [36] pipeline, with default parameters. This allowed gap filling of ambiguity regions (produced by short-read scaffolding), and further rescaffolding using long-read information. To evaluate the final assembly, several metrics and software were used. In addition to BUSCO and QUAST metrics, read back mapping of paired-end (PE) reads with Burrows-Wheeler Aligner (BWA) v.0.7.17-r1198 (BWA, RRID:SCR_010910) [37], long reads with Minimap2 v.2.17 (Minimap2, RRID:SCR_018550) [38] and RNA sequencing (RNA-Seq) with Hisat2 v.2.2.0 (HISAT2, RRID:SCR_015530) [39, 40], were also applied. To check consensus quality (QV) and k-mer completeness we used Merqury v.1.1 [41] (Figure 2).

The repetitive elements of the genome were predicted and masked by RepeatMasker v.4.0.7 (RepeatMasker, RRID:SCR_012954) [42] using homologous comparisons and ab initio predictions. First, the de novo library of repetitive elements was created with the RepeatModeler v.2.0.1 (RepeatModeler, RRID:SCR_015027) [43]. Next, the ab initio library, as well as the Dfam_consensus-20170127 (Dfam, RRID:SCR_021168) [44] and RepBase-20181026 (Repbase, RRID:SCR_021169) [45], were used in RepeatMaker to softmask the S. colias genome assembly. The genome annotation was performed with the BRAKER2 pipeline v.2.1.6 (BRAKER, RRID:SCR_018964) [46–48]. Initially, the liver RNA-Seq reads (accession number: SRR6367407 [10]) were downloaded, mapped against the S. colias genome assembly using Hisat2 v.2.2.0 [39, 40] with default parameters, and converted to BAM and sorted files using Samtools v.1.9 (SAMTOOLS, RRID:SCR_002105)  [49]. Additionally, we collected 89 proteomes from NCBI RefSeq (RefSeq, RRID:SCR_003496) [50] and Ensembl (Ensembl, RRID:SCR_002344) [51] databases. The species and accession numbers of the proteomes used in the genome annotation of S. colias can be consulted in Additional File 2 [17, 18]. Of these, 82 species belong to the class Actinopterygii (32 taxonomic orders): 81 with genome assembly at chromosome level, and 1 at scaffold level. As of the date of this genome annotation, only 1 Scombriforme genome, Thunnus orientalis, was annotated at scaffold level. The remaining 7 proteomes were selected from other vertebrate non-teleost animal models: Callorhinchus milii, Amblyraja radiata, Scyliorhinus canicula, Lepisosteus oculatus, Petromyzon marinus, Mus musculus, and Homo sapiens. Next, the RNA-Seq alignment, as well as all the above-mentioned proteomes, were inputted to the BRAKER2 pipeline with optimised parameters (–etpmode; –softmasking; –UTR = off; –crf; –cores = 30). The final file of predictions (braker.gtf) was further filtered by evidence, keeping only gene predictions with RNA-Seq or protein evidence (using BRAKER2 auxiliary scripts; selectSupportedSubsets.py), then converted to .gff3 format (using the Augustus auxiliary scripts; gtf2gff.pl) and post-processed with Another Gtf/Gff Analysis Toolkit (AGAT) v.0.6.0 [52]. The post-processing stage involved the correction of overlapping gene prediction coordinates and the removal of small or incomplete protein-coding genes (i.e., coding for <100 amino acids (aa); lacking start or stop codons). Furthermore, the proteins were extracted with AGAT and subject to functional annotation using InterProScan v.5.44.80 (InterProScan, RRID:SCR_005829) [53] and BLASTP (BLASTP, RRID:SCR_001010) searches against RefSeq  [50] and UniProtKB/SwissProt (UniProtKB, RRID:SCR_004426) [54] databases. The homology searches were performed with DIAMOND v.2.0.11.149 (DIAMOND, RRID:SCR_016071) [55] with optimised parameters (-k 1, -b 10, -e 1e-5, –ultra-sensitive, –outfmt 6). Finally, the genome and the annotation datasets were integrated using JBrowse2 (JBrowse, RRID:SCR_001004) [56], a dynamic web platform for genome visualisation and analysis that allows easy and interactive exploration of provided data (http://portugalfishomics.ciimar.up.pt/app/scombercolias/). The FASTA file containing the genome was indexed with Samtools faidx v.1.9 [49] and added to the JBrowse component, along with the annotation file sorted with GenomeTools v.1.6.1 (GenomeTools, RRID:SCR_016120) [57], and indexed with Samtools tabix v.1.9 [58]. In addition to the JBrowse component, NCBI-BLAST+ v.2.12.0 [59] was integrated into the webpage, allowing BLAST results from the genome, mRNA, protein-coding sequences (CDS), and proteins, directly from the website.

To generate a phylogenomic analysis, the proteomes of 15 selected Actinopterygii species, including the Scombriformes species Thunnus maccoyii and T. orientalis, were downloaded from public databases. The species and accession numbers used in the phylogenomic analyses can be consulted in Additional File 2 [17, 18]. Single-copy orthologs between these 15 species and S. colias were retrieved from the protein datasets by constructing protein family clusters using OrthoFinder v.2.4.0 (OrthoFinder, RRID:SCR_017118) [60] with optimised parameters (-M). This resulted in a total of 392 single-copy orthologous sequences that were individually aligned using MUSCLE v.3.8.31 (MUSCLE, RRID:SCR_011812) [61] with default parameters. Each alignment was trimmed using TrimAl v.1.2 (trimAl, RRID:SCR_017334) [62] with a gap threshold of 0.5 with optimised parameters (-gt 0.5), and afterwards concatenated using FASconCAT-G [63]. Phylogenetic inferences were conducted in IQ-Tree v.1.6.12 (IQ-Tree, RRID:SCR_017254) [64] with optimised parameters (-bb 10000 -nt AUTO -st AA). The best-fit molecular evolutionary model used in the phylogenetic analyses was JTT+F+R4, which was selected by ModelFinder [65] implemented within IQ-Tree.

To demonstrate the value of the present genome resourse, we collected the repertoire of nuclear receptors (NRs) in S. colias via TBLASTN (TBLASTN, RRID:SCR_011822) searches in the primary genome assembly with default parameters. Protein sequences of DNA-binding domains and ligand-binding domains in H. sapiens NRs were collected from the RefSeq [50] database and used in a query (NP_000466.2, NP_068804.1, NP_003241.2, XP_005257609.1, NP_001349802.1, NP_068370.1, NP_599022.1, NP_009052.4, NP_001351014.1, XP_005260464.1, NP_002948.1, NP_001257330.1, NP_003288.2, XP_016862607.1, NP_001273031.1, NP_005645.1, NP_001278159.1, NP_004442.3, NP_000167.1, XP_005268879.1, NP_004950.2, NP_201591.2). Next, regions aligning with H. sapiens sequences were collected, translated to protein using the Bio.Seq module of Biopython v.1.75 (Biopython, RRID:SCR_007173) [66], and blasted (BLASTP) against a local database containing the NR proteins of Danio rerio (D. rerio NRs database protocol: (1) NRs sequences and classifications were retrieved from [67]; (2) an NRs database was built using the makeblastdb application of NCBI-BLAST+ v.2.12.0). For each NRs sequence in S. colias, the best blast hit in the D. rerio database was collected. In some cases, several NRs of S. colias matched the same receptor in D. rerio. In these cases, the nucleotide sequences of S. colias were again validated against the NT-NCBI database, and all sequences matching different GeneIDs in the same organism were kept in the final table of NRs. In parallel, and to assess the genome annotation performed by BRAKER2, the genomic coordinates of regions aligning with H. sapiens were searched and identified in the annotation files.

To identify the genes related to the chemical defensome, target genes were selected based on a previous report profiling the “chemical defensome” of teleost species [68]. Next, gene names were used as queries to search the deduced S. colias genome annotation, a simple but successful approach for well-annotated genomes such as D. rerio [68]. When gene names were not retrieved from S. colias genome annotation (i.e., fthl, gstp, hsph, maff, nme8, slc21), further TBLASTN searches were performed in the primary genome assembly with optimised parameters (-max_hsps 1 to keep the best query-subject pair), using D. rerio sequences as a query.

To explore the variation in the demographic history of the species, a pairwise sequentially Markovian coalescent (PSMC) (PSMC, RRID:SCR_017229) strategy was applied [69], following the authors’ instructions. Briefly, PE short reads were aligned to the repeated masked genome assembly using BWA v.0.7.17-r1198 (BWA, RRID:SCR_010910) [37] with optimised parameters (BWA-MEM), and the output converted to BAM format and sorted using Samtools v.1.9 [49] (function: sort; parameters: default). Next, Picard Tools v.2.19.2 (Picard, RRID:SCR_006525) was used to remove duplicate reads (function: MarkDuplicates; parameters: default), and SAMtools for mapping quality filtering and SNP calling (function: mpileup; parameters: -Q 30 -q 30 -C 50). BCFtools v.1.9 (SAMtools/BCFtools, RRID:SCR_005227) was applied to extract consensus sequences (function: call; parameters: -c), and the subscript vcfutils (from SAMtools) was used for filtering the output for a minimum depth of 25, a maximum depth of 150, and a min RMS mapQ of 20 (function: vcf2fq; parameters: -d 25 -D 150 -Q 20). The resulting fastq file was converted to a PSMC-compatible input format using fq2psmcfa with a minimum quality threshold of 20 (parameters: -q 20). Inferences of population history were performed by running PSMC for 25 iterations with optimised parameters (-N 15, -r 5, -p 4*4 + 13*2 + 4*4 + 6) following recent PSCM estimations on Scombriformes [70]. Furthermore, to account for uncertainties in the PSMC estimates, bootstrapping of 100 replicates was performed using the split face script provided by PSMC authors. Finally, to scale the demographic estimations, a mutation rate (μ) of 7.3 × 10−9 substitutions/site/generation was used, based on a recent estimation for the Scombriformes species Thunnus albacares [70], and a generation time for S. colias of 2 years [7, 71].