Samples of Merluccius polli were collected from the waters of Mauritania (18°18’N, 16°40’W) and stored in absolute ethanol at 4°C. The morphological identification by experts was later validated by PCR amplification and Sanger sequencing of one mitochondrial (Control Region) and one nuclear (5s rDNA) marker. Later, 1µg of genomic DNA was extracted from fin tissue of a sole individual of Merluccius polli using a phenol:chloroform based protocol. The resulting DNA extraction was purified by Mag- Bind^® Total Pure NGS magnetic beads. The library was then submitted to a quality control check: DNA concentration was measured with Qubit 2.0 Fluorometer, and the integrity of the sample was checked in Bioanalyzer DNA 12000 chip. (DIN value = 8.3).

Sequencing was carried out combining Illumina NovaSeq and PacBio Sequel II sequencing in Macrogen.Inc, South Korea. SMRT libraries were prepared and subsequently sequenced with PacBio Sequel II, generating long subreads of a targeted insert size between 10 and 20 kb. PacBio subreads were assembled using Hierarchical Genome Assembly Process (HGAP 4) and default settings (Chin et al., 2013) for an estimated genome size of 534 Mb. A second library of paired-end 150 bp Illumina short reads was sequenced to correct consensus sequence generated from the PacBio data. Illumina reads with a phred score below 30 were filtered. The assemblage was corrected using Pilon v1.21 (Walker et al., 2014). K-mer analysis was employed to estimate the genome size. CDS from close species (i.e., Gadus morhua, Myripristis murdjan, Acanthopagrus latus, Archocentrus centrarchus and Anabas testudineus) were downloaded from GenBank (GCF_902167405.1, GCF_902150065.1, GCF_904848185.1, GCF_007364275.1, GCF_900324465.2) and used to make a training model using SNAP v2.31.8. Then, gene model prediction was carried out by Maker (Holt & Yandell, 2011) v2.31.8 and predictions = 1 were removed to avoid false positives. Protein prediction sets were also tested using InterProScan v5.30-69.0 (Jones et al., 2014) and psiblast v2.4.0 (Camacho et al., 2009) with EggNOG DB v4.5 (Huerta-Cepas et al., 2016). Finally, the level of completion of the genome was assessed using Benchmarking Universal Single-Copy Orthologs software (BUSCO) version 5.1.2 (Manni et al., 2021).

We conducted a bibliographic search using various search engines (i.e. Google Scholar, Web of Science) to select a set of candidate genes that, according to current scientific literature, are associated with changes in environmental factors linked to global change occurring along the distribution of Merluccius polli. Such association may be evidenced from correlational studies based on SNP, experiments for gene functions and/or other signals. The environmental factors considered were: temperature, heavy metal pollution, hypoxia, depth, and heavy metals concentration. Temperature is indeed the key factor in current climate change, and many studies describe genes involved in fish adaptation to changes in temperature (e.g. Hori et al., 2010; Dietrich et al., 2018; Lou et al., 2022). Heavy metals are main pollutants in West African waters, as explained above. The genomic mechanisms of defense of Teleostean fish against these pollutants have been described in Gadids (Olsvik et al., 2011) as well as in other species (Kim et al., 2016; Eide et al., 2021). Hypoxia is known to happen associated with the Benguela regime shift (Hutchings et al., 2009) and Benguela Niño (Monteiro et al., 2008) at the south of M. polli distribution. The genes involved in responses to hypoxia have been also described in marine fish (Tiedke et al., 2014; Xin et al., 2022). Finally, fish species change depth as a response to escape adverse conditions caused by climate change (Thresher et al., 2007; Dulvy et al., 2008), and those changes encompass a variety of physiological and genomic responses (Brown & Thatje, 2014).

A candidate gene was selected when a relation between such gene and any of the aforementioned environmental factors had been previously established on Teleosts in the literature. Genes found in phylogenetically close species (i.e. Gadiformes) were prioritized, but any Teleostean was considered a valuable example. Data were collected on the name of the gene, gene family, target factor/s to which the gene was associated, species where it was found, NCBI accession number and bibliographic reference.

To identify these candidate genes on the M. polli genome, first they were searched in the reference genes from Gadus morhua or, in case that one gene could not be found, from Danio rerio, in both cases from Ensembl version 107 (Cunningham et al., 2022). In the case that the candidate gene was identified in the source publication with a valid accession code or the proper sequence, the corresponding gene in Gadus morhua or Danio rerio was assigned by the best alignment result (the highest bit-score and lowest e-value) against the annotated protein sequences. The alignment was performed using the tool glsearch from the FASTA suite version 36 (Pearson and Lipman, 1988) taking as subject the sequence of the candidate gene using as a threshold an e-value of 1e^-15.

For the genes that had obsolete accession codes or were not found in the original publication, the corresponding gene in Gadus morhua or Danio rerio was searched by gene symbol. Next, to identify the orthologous gene in M. polli, OrthoFinder (Emms and Kelly, 2019) version 2.5. was run with default settings, specifying the parameter –f with the folder containing the protein sequences of the following species retrieved from Ensembl version 107: Danio rerio, Gadus morhua, Oreochromis niloticus, Oryzias latipes, Salmo salar and M. polli.

The candidate genes found in M. polli genome were classified by their response to environmental factors (according to literature as explained above) in four categories: responsive to temperature, depth, hypoxia, and exposure to heavy metals. Some of them may be classed in more than one category because, as part of the fish defensome, they may respond to more than one stressor (Eide et al., 2021), like regulatory transcription factors, antioxidant proteins and others.

To compute the proportion of target protein covered by its orthologous in Gadus morhua or Danio rerio, a global alignment using the Needleman-Wunsch algorithm was performed, and the proportion was defined as one minus the ratio between the number of deletions and the length of the target.

Additionally, we selected five GO functions expected to be involved in responses to the selected factors (GO:0003774, cytoskeletal motor activity; GO:0051015, actin filament binding; GO:0019825, oxygen binding; GO:0046872, metal ion binding; GO:0001666, response to hypoxia; GO:0009408, response to heat) to assess completion of our annotation in relation to functions that may potentially be linked to responses to global change. For this purpose, we searched genes associated to these GO terms in Danio rerio and Gadus morhua and matched those to our annotation of M. polli genome.

The Illumina reads generated in this study were mapped to the M. polli assembly with bwa version 0.7.17 (Li, 2013) with default parameters. The SNP calling was done using bcftools version 1.11 (Danecek et al., 2021) and default settings. The resulting SNPs were filtered so that SNPs with at least 15 reads and no more than 90 reads were considered. These values correspond to one third and twice the average coverage of the Illumina reads in the assembly.

RepeatModeler version 2.0.3 (Flynn et al., 2020) was used to build the repeat library of the M. polli assembly with the option –LTRStruct. Then, RepeatMasker version 4.1.3 was used with the Dfam database version 3.6 to scan the assembly.

This Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession JAOPHQ000000000. The version described in this paper is version JAOPHQ010000000; organism Merluccius polli C29. This new assembly was built from a combination of PacBio and Illumina sequencing data. A total of 7,349,474 subreads were obtained using PacBio Sequel sequencing accounting for 81,293,939,445 bases with a subread N50 of 16,960 (SRR21859253). From Illumina NovaSeq sequencing 157,168,555 pairs-end reads were obtained initially (SRR21859252), and after filtering the final number were 86,964,178. The details of the assembly produced are shown in Table 1 , in comparison to available assemblies to date of the genus Merluccius.

The total length of this new assembly is approximately 40% longer than previous ones ( Table 1 ). However, the most significant improvement is the level of fragmentation, which is highly reduced. The number of total sequences is down to 6,823 from more than one hundred thousand, and the statistic N50 clearly shows that the majority of the new scaffolds are longer than 160kb while in the previous assemblies, even the longer sequence was not beyond 70kb.

In order to assess the completeness of the assembly, the tool BUSCO (Manni et al., 2021) was used against the Actinopterygii lineage. The result ( Table 2 ) shows that more than 78% of the gene groups searched could be found in the assembly, with most of them complete.

The annotation carried out with Maker (Holt & Yandell, 2011) revealed a total of 26,143 genes ( Supplementary Table 1 ).

The list of candidate genes is available in Supplementary Table 2 . We identified 42, 28, 37 and 28 candidates related with the response of Teleostean fish to depth changes, exposure to heavy metals, temperature changes and hypoxia, respectively. It is worth noting that some of the candidate genes were associated with more than one factors. A total of 68 of those candidate genes were located in the new reference genome of M. polli ( Supplementary Table 2 ): 26 for depth (62% of the potential candidates), 19 for heavy metals (67.8%), 23 for temperature (62.2%) and 18 for hypoxia (64.3%).

One of the main genes involved in heavy metal responses, the metallothionein, was not found according to the annotation tool used. However, taking as a template the nucleotide sequence of this gene from Gadus morhua and using the tool blat (Kent, 2002), it is possible to identify this gene in the new M. polli in contig3709 (see the sequence in Supplementary Table 3 ). Other metal-responsive genes associated with heavy metals were indeed found, as well as genes of response to oxidative stress that are expressed in response to heavy metal pollution in other fish ( Supplementary Table 2 ).

Furthermore, the comparison of presence of genes from GO functions GO:0003774, GO:0051015, GO:0019825, GO:0046872, GO:0001666, and GO:0009408 between Danio rerio and Gadus morhua and our genome showed different percentage of representativeness among the different functions ( Supplementary Figure 1 ). GO function representativeness was overall lower for the comparison between Danio rerio and M. polli (averaging at 53.1%) than against Gadus morhua (averaging at 76.8%). This is to be expected due to the phylogenetic distance between species, as well as for the higher completion of Danio rerio genome annotation, as this species is widely used as a model organism.

From the list of candidate genes not only the majority (68 out of 109) of them can be found in this assembly, but also the predicted sequence can be considered as full-length ( Figure 2 ), taking as reference the protein sequence of G. morhua or Danio rerio used as target. The proportion of the target sequence that is covered by the identified sequence in M. polli is greater than 0.90 in 70% of them. This shows that, although the annotation was based on predictions without experimental data, the genes found are generally not fragmented.

The repeat content of this assembly is 9.24%. The amount of each type of repeats is shown in Table 3 . Of those, simple repeats are the most represented with over 7% of the genome. From interspersed repeats, that account 1% of the genome, LINEs and LTRs cover 0.5% and 0.26%, respectively. These values are much lower than related species. For example, in the Atlantic cod, gadMor2 assembly (Tørresen et al., 2017), the percentage of interspersed repeats is 23%. However, the difference is lower when looking at LINEs and LTRs where the proportion of the cod genome is 2.86% and 3.47%, respectively.

Taking advantage of the Illumina short reads generated in this study, the level of heterozygosity of the M. polli assembly was assessed. Using the standard pipeline of bcftools (Danecek et al., 2021) 680,445 heterozygous SNPs were identified, what means a heterozygosity rate of 1.16 SNPs/Kb. These rates of heterozygosity, while slightly lower than what was found in other marine fishes fall within the expected range. Other species, namely Gadus morhua (Star et al., 2011) or Larimichthys crocea (Wu et al., 2014) have higher heterozygosity rates (2.09 SNPs/Kb and 3.58 SNPs/Kb respectively). Lower heterozygosities are found in some species of sharks: both the great hammerhead and the whale shark present very low heterozygosity rates (0.53 SNPs/Kb and 0.65 SNPs/Kb respectively) (Stanhope et al., 2023). Although these heterozygosity estimates come from a single individual, this represents the first estimate of genomic variation in a species from the genus Merluccius.