Genetics in the ocean’s twilight zone: Population structure of the Mueller’s pearlside across its distribution range

María Quintela1*, Alejandro Mateos-Rivera1, Roger Lille-Langøy1, François Besnier1, Konstantinos Tsagarakis2, Naiara Rodríguez-Ezpeleta3, Miguel Bao4, Martin Wiech4, Lucilla Giulietti4, Sonia Rábade-Uberos5, Fran Saborido-Rey5, Malika Chlaida6, Espen Strand7, Sofie Knutar1, Eva García-Seoane8, Webjørn Melle7, Bjørn-Erik Axelsen7 and Kevin A. Glover1

    1Population Genetics group, Institute of Marine Research, Bergen, Norway

    2Hellenic Centre for Marine Research, Institute of Marine Biological Resources and Inland Waters, Athens, Greece

    3AZTI Basque Research and Technology Alliance (BRTA), Marine Research, Sukarrieta, Spain

    4Seafood Hazards group, Institute of Marine Research, Bergen, Norway

    5Integrated Marine Ecology, Institute of Marine Research (IIM-CSIC), Vigo, Spain

    6National Institute of Fisheries Research (INRH), 2 Avenue Sidi Abderrahmane, Casablanca, Morocco

    7Plankton group, Institute of Marine Research, Bergen, Norway

    8Sustainable Oceans and Coasts, Møreforsking AS, Ålesund, Norway

    *Corresponding author; email: maria.quintela.sanchez{at}hi.no

      bioRxiv preprint DOI: https://doi.org/10.1101/2025.02.22.639266

      Posted: February 27, 2025, Version 1

      Copyright: This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at http://creativecommons.org/licenses/by/4.0/

      ABSTRACT

      The large estimates of mesopelagic fish biomass have fuelled harvesting interests in the relatively untouched ocean’s twilight zone. The Mueller’s pearlside, one of the most abundant species inhabiting the north Atlantic mesopelagic layer, is a candidate to such fisheries despite its enormous ecological importance and the insufficient knowledge about its population genetic structure. To shed light on the latter, 863 individuals sampled across the North Atlantic and Mediterranean were genotyped using 170 genome-wide SNPs. Analyses revealed habitat-driven differentiation in three units: Mediterranean Sea, oceanic samples, and Norwegian fjords. These groups were not completely isolated to each other as a cline of Mediterranean admixture was detected in the Eastern Atlantic façade up to 47°N in an otherwise genetically homogeneous oceanic cluster. Temperature seemed to modulate the differentiation patterns, and in the Mediterranean added to the complicated topography of the Greek Seas to shape genetic structure. In the Norwegian coastline, sills did not hamper genetic exchange among fjords ranging 200 km apart, probably due to the position of the species in the water column together with its swimming capacity. This genetic information should be combined with demographic properties to outline the management of this species prior to any eventual fishery attempt.

      INTRODUCTION

      The mesopelagic zone, located between 200–1000 m depth, and hosting up to 90% of the total oceanic fish biomass (Irigoien et al., 2014), is a largely unexplored area. The limited amount of sunlight at these depths explains the synonym ocean’s twilight zone, driving unique adaptations such as bioluminescence and diel vertical migration in the organisms inhabiting it. The enormous biomass of mesopelagic fish, ambitiously estimated around 10 Gt (10,000 million tonnes) (Irigoien et al., 2014), fuelled interests for the exploitation of a resource intended to aid meeting the demands of the expanding world population, both in terms of human food security and animal feed in a time where over a third of the worldwide fisheries operate beyond their biologically sustainable levels (FAO, 2020). However, a sustainable exploitation of the mesopelagic resources encounters challenges such as the limited knowledge about trophic interactions, life histories, behaviour, biomass, diversity, nutritional composition, and population genetic structure of the mesopelagic organisms (Hidalgo & Browman, 2019Alvheim et al., 2020Martin et al., 2020Standal & Grimaldo, 2020Wiech et al., 2020).

      Mesopelagic organisms play a key ecological role in the Biological Carbon Pump (BCP), which is the suite of processes that collectively account for the ocean’s biologically driven sequestration of carbon from the atmosphere and land runoff to the ocean interior and seafloor sediments (Ducklow, 2001Sigman & Haug, 2006Boyd et al., 2019). Photosynthesis by phytoplankton in the sunlit surface ocean involves net uptake of CO2 from the atmosphere, whereas the transference of carbon from the surface to the deep ocean, amounting 1,300 Pg C yr-1, occurs via distinct pathways including gravitational settling of organic particles, mixing and advection of suspended organic carbon, as well as active transport by vertically migrating metazoans (Nowicki et al., 2022). The diel vertical migrations (DVM) of many mesopelagic fish species that feed on zooplankton near the surface at night and return to deep layers during the day actively export carbon from the surface to deeper water masses and contribute to the biological pump (Hidaka et al., 2001Shreeve et al., 2009Robinson et al., 2010) in what is regarded as the “largest daily migration of animals on earth” (Hays, 2003). The active transport mediated by diel vertical migration sequesters more carbon than the physical pump (subduction + vertical mixing of particles), i.e. 1.0 vs. 0.8 Pg C because of deeper remineralization depths (Stukel et al., 2023). Without it, the levels of atmospheric CO2 would be about 400 ppm higher than at present (Sanders et al., 2014Boyd, 2015).

      The Mueller’s pearlside, Maurolicus muelleri (Gmelin, 1789) is a small but abundant mesopelagic fish belonging to the family Sternoptychidae (see Grimaldo et al. (2020)). The genus Maurolicus was initially regarded as monotypic with a single species, M. muelleri (Gmelin, 1789), of cosmopolitan distribution. Later on, following identification of differences in combinations of morphometric characters, the genus was split into 15 species (Rees et al., 2020, see Fig. 1). However, many of them display overlapping ranges of meristic and morphometric traits and thus remain poorly characterized (Rees et al., 2017). Mitochondrial (16S and COI) and nuclear (ITS-2) gene sequences revealed groups conflicting with previously recognised species: (1) a ‘Northern’ clade comprising M. muelleri and M. amethystinopunctatus, (2) a ‘Southern’ clade comprising M. australisM. walvisensis (also M. japonicus) and (3) an Eastern Equatorial and Western North Atlantic M. weitzmani (Rees et al., 2017). Synonymisation is proposed for M. muelleri and M. amethystinopunctatus, with limited morphological variation likely to reflect physical and biological differences experienced North/South of the sub-polar front (Rees et al., 2017). In a follow-up involving multiple locations worldwide, M. muelleri and M. australis were regarded as potentially a single species, but since no shared haplotypes were found between the two disjunct groups, they were kept separately (Rees et al., 2020). M. muelleri is thus distributed in the North Atlantic and Mediterranean Sea.

      Fig 1.

      Maurolicus muelleri samples. Blue dots depict the Ocean samples whereas green dots depict the Mediterranean ones. The fjord samples, zoomed in the box to the right, are signalled with red dots. The main cold ocean currents are outlined in blue whereas the warm ones are outlined in red. For three of the sampling sites (i.e. Iceland, Cantabrian Sea and Gulf of Biscay) consisting of different nearby stations, the geographic coordinates correspond to the most central site.

      M. muelleri is one of those species conducting DVM that play a critical role in marine ecosystems both as part of the BCP as well as linking primary consumers and higher tropic levels such as larger fish, cephalopods, seabirds and mammals (e.g. Cherel et al., 2010Drazen & Sutton, 2017). The species forms distinct sound scattering layers at depth (SLs) and can be acoustically detected even from the youngest ontogenic stages due to their prominent swim bladder. In the Norwegian fjords, the SL composed of adults display DVM between April and June (Rasmussen & Giske, 1994Goodson et al., 1995), whereas in early winter (January) adult fish remain in deeper waters while juveniles perform DVM as well as nocturnal midnight sinking after dusk (Giske et al., 1990Baliño, 1993Goodson et al., 1995). Fully grown individuals can achieve swimming speeds of 2–7 body lengths per second (i.e. 8–30 cm·s-1) (Torgersen & Kaartvedt, 2001).

      In addition to taxonomic considerations, the literature dealing with population structure in M. muelleri is sparse and limited to few genetic studies. In western Norway, allozymes revealed that samples collected from five fjords displayed varying degrees of genetic differentiation to each other, and to a sample collected in the North Sea (Suneetha & Nævdal, 2001). In the Bay of Biscay, thousands of SNPs obtained from RADseq revealed no clear evidence of population structure as well as lower observed heterozygosity than anticipated (Rodriguez-Ezpeleta et al., 2017). Finally, mitochondrial DNA genes (COI, 12S, and 16S) revealed a lack of genetic structure among samples taken from the Greek Seas (Sarropoulou et al., 2022).

      The correct outline of biologically correct management units or stocks is one of the multiple requirements of sustainable fisheries and intends to prevent the overexploitation of unique spawning components (e.g., see Kerr et al. (2017) for review). As this ecologically crucial species is a potential candidate for exploitation in a commercial fishery, there is a need to fill the knowledge gap regarding its ocean-wise genetic structure. In addition, when outlining management areas, special attention needs to be directed to identify locally adapted populations or subdivided ones as they might have different sustainable yield levels and be more prone to the negative effects of overfishing (Waples, Punt, and Cope 2008; Pinsky and Palumbi 2014). Furthermore, the range shifts that many marine species are already experiencing in the face of climate change should also be taken into consideration (Palacios-Abrantes et al. 2022; Dahms and Killen 2023).

      The aim of the present study was to provide the first ocean-wide examination of genetic structure of M. muelleri. Several hundreds of individuals sampled across most of the species’ distribution range including both sides of the Atlantic (between longitudes 26 °W–25°E and latitudes 35–70 °N) were high-throughput genotyped using a suite of 170 SNP loci. Using the ecological parallelism with the mesopelagic glacier lanternfish Benthosema glaciale (Quintela et al., 2024), we hypothesize that the genetic patterns of differentiation of M. muelleri could be habitat-related with three major units corresponding to Norwegian fjords, open Atlantic ocean and Mediterranean, respectively.

      MATERIALS AND METHODS

      Sampling and genotyping

      In the period 2017–2024, 941 individuals were collected in 23 sampling sites ranging across large part of the North Atlantic distribution of the species (Fig. 1). Based upon our previous experience with another widely distributed mesopelagic fish, B. glaciale (Quintela et al., 2024), three main habitats were aimed at: Norwegian fjords (i.e. five samples around 60 °N), Mediterranean Sea (i.e. one sample from the Alborán Sea plus four Eastern samples from the Greek Seas), and open sites in the Atlantic ocean (i.e. thirteen samples within latitudes 35 – 70 °N and longitudes 26 °W–25 °E). Wherever available, 50 individuals were collected per location. Two of the Greek samples came from the enclosed, deep and isolated gulfs of Corinth (Ionian Sea) and of North Euboea, with geophysical features similar to the Norwegian fjords (Kapelonis et al., 2023). Fin clips were taken and stored in ethanol 96% prior to DNA isolation, which was conducted using SPRI paramagnetic beads from the Beckman Coulter DNAdvance kit (A48706). DNA concentration was quantified using NanoDrop 8000.

      Four sampling sites belonging to the different habitats were selected for SNP mining following the procedure used in Quintela et al. (2024): Mediterranean (Ionian Sea), fjords (Boknafjord) and oceanic (two sites in the Atlantic Ocean at latitudes 40 and 50 °N, i.e. Atlantic_40N and Celtic Sea, respectively). From these samples RNA-free DNA was extracted using the Beckman Coulter DNAdvance kit (A48706). DNA integrity was assessed by agarose gel electrophoresis and DNA concentration was quantified using Thermo Fisher Qubit dsDNA Broad Range (Q32853). Equimolar amounts of DNA from 10 individuals per site were pooled, and one library was prepared per pool using Illumina DNA prep. Pooled samples were then sequenced on a NovaSeq X using 1/8 of a NovaSeq X Series 25B flow cell (150 PE). The package fastp (Chen et al., 2018) was used for data preprocessing, which included deduplication, adaptor removal and analysis of over-represented sequences. FASTQC v0.11.5 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and MULTIQC v1.7 (https://github.com/MultiQC/MultiQC) were used to output graphics and statistics of the data quality control.

      A draft assembly was produced using MEGAHIT (Li et al., 2015Li et al., 2016) testing Kmer from 41 to 73. Scaffolding was performed with SoapDeNovo (https://github.com/aquaskyline/SOAPdenovo2) using 63 as Kmer value. The obtained draft assembly consisted mostly of small contigs, shorter than 500bp. However, 2156 contigs were larger than 10kb and represented together 27Mb. The four paired-end pooled samples were then mapped against this 27Mb short genome using BWA V0.7.17 (Li & Durbin, 2010) with default parameters for BWA aln and BWA sampe. Variant sites were called with the mpileup function (Li, 2011) from samtools V1.9 (Li et al., 2009). The R package vcfR (Knaus & Grünwald, 2016Knaus & Grünwald, 2017) was used to visualise the distribution of coverage and mapping quality across all variant sites. Only Single Nucleotide Polymorphic (SNP) variants were retained if they fulfilled the following criteria: phred-scaled quality score (QUAL >500), coverage larger than 120X and lower than 300X for each pooled sample sites, and both alternative and reference allele represented in at least two samples out of four. This filtering produced 12,000 SNP that revealed significant differentiation among the three putative habitats as hypothesized. The differentiation was particularly large versus the Mediterranean according to the first axis of the PCA plot (PC1), which accounted for 39.5% of the variation (Suppl. Fig. S1) whereas PC2 (31.2%) depicted the differentiation between fjord and ocean as well as between oceanic samples. Finally, SNPs located less than 200 bases from another polymorphic site (SNP or indel) were removed to keep only SNPs with stable primer sequences. A total of 2034 SNPs fulfilled those criteria; however, the goal was to genotype as many individuals as possible across the species’ distribution range using a high-throughput SNP genotyping approach in a reasonably affordable manner. Therefore, the Assay Designer option of the Typer v5 program (Agena Biosciences CA, USA), was used to arrange the 2034 SNPs into different multiplex reactions with a maximum of 25 loci each. Based on number of SNPs and the multiplex confidence score, eight multiplexes containing in total 200 SNPs were finally selected to genotype the full set of samples (N=941). Primers were designed, and genotyped on the 941 individuals using the Sequenom MassARRAY iPLEX Platform as described by Gabriel et al. (2009). From the Flemish Cap (NW Atlantic), only 10 individuals were obtained, which were sampled both for fin clips and muscle. Thus, DNA was extracted from both tissues and individuals were amplified twice to conduct positive controls and reproducibility.

      Genetic structure

      Statistical analyses were restricted to the subset of 170 well-functioning polymorphic loci obtained out of the 200 screened ones (Suppl. Table S1 compiles the reasons for discarding the 30 ones). Even though the threshold of acceptance of missing data per individual was 20%, only 3 of the 863 finally retained individuals showed the maximum missing data allowed, whereas 63% of them displayed ≤5%. These 863 individuals were distributed into 23 samples across three different habitats: fjords, Atlantic Ocean and Mediterranean Sea. To assess if the 170 SNPs would accurately discriminate between individuals in a population, the genotype accumulation curve was built using the function genotype_curve in the R (Team, 2020) package poppr (Kamvar et al., 2014) by randomly sampling x loci without replacement and counting the number of observed multilocus genotypes (MLGs). This repeated r times for 1 locus up to n-1 loci, creating n-1 distributions of observed MLGs. The observed (Ho) and unbiased expected heterozygosity (uHe) as well as the inbreeding coefficient (FIS) were computed for each sample with GenAlEx v6.1 (Peakall & Smouse, 2006). Likewise, the genotype frequency of each locus and its direction (heterozygote deficit or excess) was compared with Hardy-Weinberg expectations (HWE) using the program GENEPOP 4.0.6 (Rousset, 2008) as was linkage disequilibrium (LD) between pairwise loci. The False Discovery Rate (FDR) correction of Benjamini and Hochberg (1995) was applied to p-values to control for Type I errors. Non-parametric Kruskal-Wallis’s rank sum test followed by post-hoc Dunn’s test were applied to perform comparison of estimates of genetic diversity among habitats. Data conversion into the appropriate formats was conducted using PGDSpider 2.1.1.5 (Lischer & Excoffier, 2012).

      Unsupervised analysis of genetic structure were conducted with Principal Component Analysis (PCA) using the function dudi.pca in ade4 package (Dray & Dufour, 2007) in R (Team, 2020) after replacing missing data with the mean allele frequencies, using no scaled allele frequencies (scale = FALSE). In addition, the Bayesian clustering approach implemented in STRUCTURE v.2.3.4 (Pritchard et al., 2000), and conducted using the software ParallelStructure (Besnier & Glover, 2013), was used to identify genetic groups under a model assuming admixture and correlated allele frequencies without using LOCPRIORS. Ten runs with a burn-in period consisting of 100,000 replications and a run length of 1,000,000 MCMC iterations were performed for K=1 to K=10 clusters. To determine the number of genetic groups, STRUCTURE output was analysed using two approaches: a) the ad hoc summary statistic ΔK of Evanno et al. (2005), and b) the Puechmaille (2016) four statistics (MedMedK, MedMeanK, MaxMedK and MaxMeanK), both implemented in StructureSelector (Li & Liu, 2018). Finally, the ten runs for the selected Ks were averaged with CLUMPP v.1.1.1 (Jakobsson & Rosenberg, 2007) using the FullSearch algorithm and the G’ pairwise matrix similarity statistic, and graphically displayed using bar plots.

      Supervised genetic structure using geographically explicit samples was assessed using the Analysis of Molecular Variance (AMOVA) and pairwise FST (Weir & Cockerham, 1984), both computed with Arlequin v.3.5.1.2 (Excoffier et al., 2005). Furthermore, the relationship among samples was examined using the Discriminant Analysis of Principal Components (DAPC) (Jombart et al., 2010) implemented in the R (Team, 2020) package adegenet (Jombart, 2008) in which groups were defined using geographically explicit locations. To avoid overfitting, both the optimal number of principal components and discriminant functions to be retained were determined using the cross-validation function (Jombart & Collins, 2015Miller et al., 2020).

      The relationship between genetic (FST) and geographic distance was examined to investigate if it adhered to the expectations of an “Isolation by Distance” pattern (IBD), i.e. increasing genetic differentiation with geographic distance as a result of restricted gene flow and drift (Wright, 1943Slatkin, 1993Rousset, 1997). A two-tailed Mantel (1967) test was conducted using PASSaGE v2 (Rosenberg & Anderson, 2011) and significance was assessed via 10,000 permutations. The matrix of pairwise shortest geographic distance by water (i.e. avoiding land masses) was calculated with the R (Team, 2020) package marmap (Pante & Simon-Bouhet, 2013).

      Putative clines of allele frequency extending from the easternmost Atlantic coast into the Mediterranean were investigated via the latitudinal sliding-window approach developed by Pereira et al. (2018) using the R-scripts provided by the authors in Table S1.3. In a second step, loci identified as displaying clines were subjected to a geographic cline analysis conducted using the R package HZAR (Derryberry et al. 2014) over a circa 8000 km transect starting in Vesterålen and finishing in the Cretan Sea. The fifteen models implemented in HZAR were fitted to the allele frequency of the candidate loci to determine the position, width and shape of cline over the total geographic distance. The reference cline was built using STRUCTURE Q-score and, in both cases, the best cline model was decided upon AIC scores.

      Outlier detection and environmental association analyses

      Loci departing from neutrality were identified using a combination of Arlequin v.3.5.1.2 (Excoffier et al., 2005) and BayeScan 2.1 (Foll & Gaggiotti, 2008). Analyses on Arlequin were conducted based on 100 demes with 50,000 simulations under a hierarchal island model. In BayeScan, sample size was set to 10,000 and the thinning interval to 50. Loci with a posterior probability over 0.99, corresponding to a Bayes Factor >2 (i.e., “decisive selection” (Foll & Gaggiotti, 2006), were retained as outliers. The outcome of both methods was intersected to identify a panel of eventual consensus outliers.

      Adaptation to local environments often occurs through natural selection acting on a large number of loci, each having a weak phenotypic effect. Temperature, salinity and dissolved oxygen are variables of physiological importance for marine ectotherms and of significance to account for myctophid distributions (e.g., Koubbi et al., 2011Flynn & Marshall, 2013). Data on summer temperature measured at 200 m depth using a grid of ¼ ° and averaged for the period 2005-2012 was retrieved from NOAA database (National Oceanic and Atmospheric Administration). LFMM, “latent factor mixed model” (Frichot et al., 2013), was used to assess which variables could be a potential selective pressure driving local adaptation by identifying loci showing unusual associations with these environmental factors compared to the genetic background. This method uses a linear mixed model to test for associations between genetic variation and environmental factors, while controlling for neutral genetic structure with (random) latent factors. The number of latent factors was set at K=3 according to DAPC outcome and ten runs were conducted using 1,000 sweeps for burn-in and 10,000 additional sweeps. The corresponding z-scores of the ten replicates were combined following the recommendations described in Frichot and François (2015). First, the genomic inflation factor (λ) was obtained after computing the median of the squared (combined) z-scores for each K, divided by the median of the χ2 distribution with one degree of freedom. Finally, p-values were adjusted using the genomic inflation factor (λ) and false discovery rates were set using the Benjamini and Hochberg (1995) algorithm.

      Redundancy Analysis (RDA), a genotype-environment association (GEA) method to detect loci under selection (Forester et al., 2018), was be used to analyse loci and several environmental predictors simultaneously. The SNP genotype matrix was used as the response matrix, while the environmental data matrix served as the explanatory variables. RDA determines how groups of loci covary in response to the multivariate environment, and can detect processes that result in weak, multilocus molecular signatures (Rellstab et al., 2015Forester et al., 2018). Environmental data (temperature, salinity, pH, dissolved oxygen, current velocity, chlorophyll) was obtained for the different sampling points using the Bio-ORACLE database https://www.bio-oracle.org (Tyberghein et al., 2012Assis et al., 2018). Collinearity between variable pairs was investigated and only non-correlated ones were retained for analyses. Analysis was conducted with the R package vegan v.2.5–7 (Oksanen et al., 2020). The candidate outliers obtained from the different analyses were annotated by matching the SNP flanking regions against non-redundant database of GenBank (www.ncbi.nlm.nih.gov/genbank/) using the Basic Local Alignment Search Tool (Altschul et al., 1990).

      RESULTS

      Genetic structure

      SNPs genotyping of the 10 fish collected in the Flemish Cap which DNA was extracted in duplicate revealed similar performance for muscle and finclips, except in one of the individuals where finclips failed in 92% of the 200 tested loci. In the remaining nine pairs and excluding the differences due to failure of amplification in one of the individuals, 14 discrepancies (i.e. homozygote vs. heterozygote) were found (0.77%).

      When considering the full set of individuals, the plateau of the genotype accumulation curve reached with 17 out of the 170 retained loci implied that 10% of the SNPs displayed enough power to discriminate between unique genotypes from the same sample (Suppl. Fig. S2). The Kruskal-Wallis test detected significant differences in genetic diversity across habitats (P of 0.005 and 0.002, respectively for HO and uHE, and P=0.05 for the percentage of polymorphic markers) with lower values reported for the Mediterranean in all cases (Table 1). Dunn’s a posterior test revealed significant differences for the three estimates between Mediterranean and Oceanic habitats and between Fjords and Mediterranean for HO. Out of the 330395 performed tests for LD, 3% were significant after FDR. Deviations from HWE were found in 49 out of the 3910 loci by population tests after FDR (1.25%).

      Table 1.

      Sample summary statistics obtained for the set of 170 SNP loci: Sampling sites arranged into habitats with geographic coordinates in decimal degrees; number of individuals (N); percentage of polymorphic loci (Pol. loci); observed heterozygosity, Ho (mean ± SE); unbiased expected heterozygosity, uHe (mean ± SE); inbreeding coefficient, FIS (mean ± SE); number of deviations from Hardy-Weinberg equilibrium (Dev. HWE) and number of deviations from Linkage Disequilibrium (Dev. LD) at α=0.05 both before and (after) False Discovery Rate (FDR) correction. As some of the sampling sites consists of different nearby stations, the geographic coordinates indicated per sample are an average of all trawls for simplicity.

      The first axis of the PCA biplot, accounting for 20% of the variation, neatly separated the Mediterranean from the fjords and part of the Oceanic individuals; however, a transition could be observed between the N. Euboean Gulf and some of the Atlantic samples such as Morocco, the Cantabrian Sea or the Gulf of Biscay (Fig. 2). The second axis (3.8%) separated the fjords and the Oceanic samples. In STRUCTURE a posteriori analyses, the Evanno test strongly suggested K=2 as the most likely number of genetic groups (ΔK=33400, Suppl. Fig. S3). This first division singled out the Mediterranean samples while revealing a decreasing cline of the Mediterranean genetic profile into the oceanic samples from Morocco northwards the Gulf of Biscay (Fig. 3a). The uniformity of the Mediterranean samples was disrupted in the Greek sample of the N. Euboean Gulf. At K=3, fjord and oceanic clusters became evident while the cline Mediterranean-Atlantic was still present in the southern Atlantic samples (Suppl. Fig. S4a). The Mediterranean sample from N. Euboean Gulf revealed a different cluster at K=5 (Suppl. Fig. S4c). Puechmaille’s solution of K=6 (Fig. 3b) reflected the genetic uniformity of the fjord samples, the distinctness and rather uniform Mediterranean genetic profile (excluding the N. Euboean Gulf) as well as two admixed clusters in the Atlantic Ocean (a northern one from Norway 63 °N to Flemish Cap and a southernmost one ranging from 47 °N to Morocco). Interestingly, from K=3 onwards, the sample from Vesterålen deviated from the oceanic profile and displayed a seemingly physical mixture of fish belonging to the fjord (majority) and northern Atlantic clusters.

      Fig 2.

      Principal Component Analysis (PCA) of Maurolicus muelleri genotyped at 170 loci. Blue dots depict the Ocean samples, red dots indicate fjord individuals and green dots represent the Mediterranean ones.

      Fig 3.

      Barplot representing the proportion of individuals’ ancestry to cluster at a) K2 and b) K7 after Bayesian clustering in STRUCTURE as determined by the a posteriori analyses conducted using Evanno’s test and Puechmaille’s statistics, respectively.

      High levels of genetic differentiation were detected on a per-locus basis, with 150 out of 170 loci significantly different from zero after FDR, and 83 of them with FST values ranging between 0.1 and 0.8. Hierarchical AMOVA revealed significant structure both among habitats (FCT=0.158, P<0.001), among samples within habitats (FSC=0.057, P<0.001) and within samples (FST=0.206, P<0.001) with 15.8% of the variation hosted among habitats and 79.4% within samples. The dendrogram coupled with the pairwise FST (Fig. 4) revealed a first dichotomic division separating the Mediterranean, Morocco and the Cantabrian Sea from all the remaining sites. Interestingly, not only did the sample from the N. Euboean Gulf cluster with the latter Atlantic samples, but it also showed lower genetic differentiation towards both the Atlantic and the fjords than any of the Mediterranean samples (Suppl. Table S2). In the second branch of the tree, the fjords formed their own ramification while joining to Vesterålen whereas the Ocean samples were divided into the northernmost plus westernmost ones versus the ones between 47 and 40 °N. Most of the pairwise FST values were significantly different from zero (Suppl. Table S2) and strong differentiation was found among habitats, particularly towards the Mediterranean but exceptions were detected within habitats. The only significant differentiation found among the fjords occurred between the two farthest most sites, i.e. Osterfjord and Boknafjord separated by 200 km (FST=0.007, P=0.001). In the Mediterranean, the only samples that did not differ from each other were the N. Aegean Sea and the Cretan Sea. The complicated topography of the Greek Seas with very confined body waters such as in the N. Euboean Gulf and the Gulf of Corinth (Ionian Sea) seems to hamper the genetic exchange whereas the 462 km of distance between the N. Aegean Sea and the Cretan Sea sampling sites did not seem to be a limitation for the gene flow. Within the oceanic samples, 71% of the pairwise comparisons were significant. However, almost no differentiation was detected between the samples ranging from Vesterålen southwards to the Celtic Sea together with the three westernmost ones (Iceland, AtlanticXIIc and FlemishCap). In the group of oceanic samples ranging between Morocco northwards to 47 °N, both the Cantabrian Sea and Morocco were significantly different from the remaining (pairwise FST ranging between 0.005 and 0.025) whereas no differentiation was detected among Atlantic_40 °N, Atlantic_47 °N and the Gulf of Biscay.

      Fig 4.

      Heatmap FST coupled with dendrogram. Pairwise values and significance can be found in Suppl. Table S2. Colours depict the degree of differentiation: beige colours indicating lower differentiation and moving towards dark brown to indicate larger differentiation.

      The DAPC revealed both differences among the three habitats (Ocean, Fjords and Mediterranean) as well as a certain continuity among them. Thus, the first axis of differentiation (60.8% of the variation) showed the uniqueness of the fjord individuals as well as a continuum between the Atlantic samples and the Mediterranean, together with the total lack of overlapping between the Mediterranean and the Northernmost and Westernmost Atlantic ones (Fig. 5a). The second axis of differentiation (16.8%) further singled out the fjord individuals. As formerly seen, some of the individuals from Vesterålen clustered with the fjord samples whereas the remaining clustered with the Northern group of the Atlantic. The third axis of differentiation (7%) singled out the sample from the N. Euboean Gulf (Fig. 5b). The positive and strong correlation (p=0.783, P<0.001) between the shortest water distance between samples and the genetic distance measured as pairwise FST was indicative of a pattern of Isolation by Distance (Fig. 6).

      Fig 5.

      Genetic differentiation among Maurolicus muelleri samples assessed with 170 SNP loci using Discriminant Analysis of Principal Components (DAPC) after retaining 100 principal components and 3 discriminant functions: a) axis 1 and 2, and b) axis 1 and 3. Individuals from different sampling sites are represented by coloured dots.

      Fig 6.

      Mantel test showing significant correlation between shortest geographic distance by water (km) and genetic distance (FST). Colours depict the samples involved in the comparisons within (red, blue, green) and across habitats (grey).

      The allele frequencies of 72 out of the 170 loci showed signs of a latitudinal cline and only in one of them, the difference between maximum and minimum frequencies was <0.2 (in 49 of the loci, delta frequency ranged between 0.5 and 1). STRUCTURE Q-score also adhered to a clinal pattern with delta frequency of 0.99 (see Suppl. Fig. S5). HZAR revealed that for 18 of the 72 loci, the cline was centred within the limits of the reference cline, i.e. around 3000–3400 km south from Vesterålen (Fig. 7), which corresponds to a latitude around 44–47 °N.

      Fig 7.

      Summary of HZAR analyses: Cline centre and width for the loci showing signs of latitudinal clines. Red dotted lines represent the limits of the reference cline built with STRUCTURE Q-score.

      Outlier detection and environmental association analyses

      Arlequin and BayeScan jointly flagged 20 loci as candidates to positive selection, 25 as putatively under balancing selection and 78 that did not depart from neutral expectations. The remaining 47 loci were found to be under selection only by one of the procedures. Data were divided into putatively identified neutral loci and loci under positive selection, and genetic structure was assessed separately.

      In the neutral dataset, both the PCA biplot (Suppl. Fig. S6a) and the DAPC (Suppl. Fig. S6b) depicted the differentiation among habitats revealed by the set of total markers but with reduced capacity of discrimination of the transition between the Atlantic and the Mediterranean. Evanno test again suggested K=2 as the most likely number of genetic groups (ΔK=2105, Suppl. Fig. S7) and the resulting barplot mirrored the result obtained with the full set of loci (Suppl. Fig. S8a). The Puechmaille’s solution of K=6 (Suppl. Fig. S8b) revealed a similar pattern to the one with total loci without discriminating the Mediterranean sample of the N. Euboean Gulf. Hierarchical AMOVA revealed significant structure both among habitats (FCT=0.127, P<0.001), among samples within habitats (FSC=0.030, P<0.001) and within samples (FST=0.154, P<0.001). The 12.7% of the variation was hosted among habitats whereas 84.6% was hosted within samples. Pairwise FST (Suppl. Fig. S9, Table S3) followed the patterns of the total loci.

      The set of 20 candidate loci to positive selection lacked capacity of discriminating between fjord samples and oceanic ones, namely the Northern Atlantic cluster, both at PCA (Suppl. Fig. S10a) and the DAPC (Suppl. Fig. S10b). This lack of discrimination is illustrated by the three main groups of the FST dendrogram (Suppl. Fig. S11, Table S4): i) Norwegian fjords and Atlantic samples north of 47°, ii) Atlantic samples from 47 °N southwards plus the N_EuboeanGulf, iii) rest of Mediterranean samples.

      LFMM analysis after adjusting p-values using the genomic inflation factor (λ) revealed that 21 loci were associated with summer temperature at 200 m depth. Of those loci, 19 were flagged as candidate outliers to positive selection both by BayeScan and Arlequin whereas the two remaining were only flagged by BayeScan. All loci identified by LFMM were involved in allele frequency clines (see Table 2 and Suppl. Table S5) with a break point at 47°N in the eastern part of the Atlantic. Interestingly, the allele frequency at locus P402 and P676 in N. Euboean Gulf fit better the frequency of the fjord and northern Atlantic samples than the Mediterranean ones, pattern that happens in some other loci and accounts for the distinctness of the N. Euboean sample and its relatedness to the Atlantic ones (see Fig. 3 or Fig. 5a). RDA identified two loci putatively related with temperature (P942 and P1458), both candidate outliers to positive selection, identified by LFMM as linked to temperature and with clines centered around 3230–3251 km. In addition, RDA identified a locus seemingly related to chlorophyll concentration (P452), which was flagged as under positive selection by BYSC and was not involved in any cline (Suppl. Table S5). BLAST identified two matches into zinc finger proteins, two in metalloproteinases, anoctamin and MAX interactor, respectively (Suppl. Table S5).

      Table 2.

      LFMM analysis: Heatmap of allele frequency per sample for the candidate loci associated to summer average temperature measured from 2005 to 2012 at 200 m depth. Samples are distributed in blocks corresponding to habitats, i.e. fjords, oceanic and Mediterranean.

      DISCUSSION

      The present study fills the knowledge gap on the population genetic structure of a prospective fisheries target, the Mueller’s pearlside Maurolicus muelleri, across most of its distribution range. Three main genetic groups inhabiting the Mediterranean Sea, the Norwegian Fjords and the Atlantic Ocean, respectively supports our hypothesis of three habitat-driven units.

      Genetic markers

      The genetic markers used in the current study stem from a two-step procedure: first, mining SNPs from pool sequencing data, and then second: genotyping a subset of selected loci in several hundred individuals using a relatively affordable and high-throughput SNP genotyping method. This same methodological approach was recently used to mine 121 SNPs that revealed remarkable differentiation across habitats in the glacier lanternfish B. glaciale (Quintela et al., 2024). In a similar manner, but using ddRAD sequencing, 91 SNPs mined for European sprat unravelled patterns of differentiation (Quintela et al., 2020) later confirmed by full genome sequencing (Petterson et al., 2024). The success of former experiences, the confirmation of the genetic patterns observed with full genome sequencing with the set of 170 selected SNPs together with the outcome of the Genotype Accumulation Curve (Supplement Fig. S2), made us confidently rely on the discriminating capacity of the SNP array used here.

      Patterns of genetic differentiation

      The genetic patterns of differentiation of M. muelleri align to a great extent with the ones observed in other mesopelagic fish of similar size, B. glaciale (Quintela et al., 2024), and add to the body of literature of genetically-identified ecotypic differentiation detected in marine taxa and driven by environmental factors such as bathymetry as in tusk (Brosme brosme) (Knutsen et al., 2009) or beaked redfish (Sebastes mentella) (Benestan et al., 2021), salinity as in European sprat (Sprattus sprattus) (Quintela et al., 2020Pettersson et al., 2024), or the distinction between marine vs. coastal habitat as demonstrated in European anchovy (Engraulis encrasicolus) (Le Moan et al., 2016), long-snouted seahorse (Hyppocampus guttulatus) (Riquet et al., 2019Meyer et al., 2024), northern shrimp (Pandalus borealis) (Hansen et al., 2021) or Atlantic cod (Gadus morhua) (Knutsen et al., 2018).

      A major difference between M. muelleri and B. glaciale is that in the latter, loci in linkage disequilibrium drive a three clusters’ pattern detectable through PCA that matches the genetic signature associated to chromosomal inversions. The arrangement of this putative inversion shows remarkable frequency differences between the Atlantic and the Mediterranean hindering the recombination after gene flow and thus preventing the development of an admixture gradient in that part of the genome. No signature of Mediterranean genetic profile is detected in the genetically uniform Atlantic Ocean samples. However, while we cannot rule out the existence of small chromosome inversions in the genome of M. muelleri that could have eluded the scrutiny of our markers, no striations were detected in the PCA and no patterns of differentiation corresponding to different haplotypes were revealed through unsupervised STRUCTURE analyses either. The lack of an apparent inversion in M. muelleri is manifested in the Eastern Atlantic samples through an allele frequency cline extending from North Moroccan waters northwards 47 °N. Clinal patterns of differentiation between habitats are also known to the transition Baltic Sea-Atlantic ocean (reviewed in Johannesson & Andre, 2006).

      In connection with this putative inversion, the Mediterranean samples of B. glaciale displayed less than half the genetic variation compared to the samples from the Atlantic and Norwegian fjords, whereas in M. muelleri the Mediterranean samples displayed slightly lower genetic diversity than the rest. This finding highlights the long-term isolation of the Mediterranean unit, which may have led to an adaptation process that reduced their genetic composition to include only alleles suitable to the environment. The values of observed and expected heterozygosity found in this study were similar to B. glaciale (Quintela et al., 2024). In contrast, Rodriguez-Ezpeleta et al. (2017) reported exceptionally low levels of genetic diversity (some 4-fold lower) in M. muelleri from the Bay of Biscay genotyped at thousands of RAD-seq SNPs, summoning the Wahlund effect. In our study, an extremely low but significant differentiation was detected between the Cantabrian Sea and the Bay of Biscay, while similar patterns were seen in STRUCTURE barplots consistent with the findings of Rodriguez-Ezpeleta et al. (2017). The geographic range of our study allows us to suggest that at K=2, the influence of the Mediterranean genetic profile could account not only for the substructure proposed by Rodriguez-Ezpeleta et al. (2017), but would seem to drive the patterns of differentiation in the Atlantic cluster in a similar manner as the introgression with North East Artic Cod shapes population structure in Norwegian Coastal Cod (Dahle et al., 2018Jorde et al., 2021). Genetic breaks between the Atlantic and the Mediterranean have been identified in multiple taxa of fish (Bargelloni et al., 2003Tine et al., 2014Barros-García et al., 2020Quintela et al., 2020) including seahorses (Riquet et al., 2019Meyer et al., 2024), as well as in sponges (Riesgo et al., 2019), molluscs (Pérez-Losada et al., 2002Lapègue et al., 2023), crustaceans (Reuschel et al., 2010Jenkins et al., 2019) and harbour porpoise (Fontaine, 2016).

      In M. muelleri oceanic samples, the Northern Atlantic cluster displays almost no differentiation across the Atlantic from Norway 63 °N to Flemish Cap (3570 km distance in a longitude range from 3.78 °E to 45 °W), which closely aligns with the lack of structure displayed by B. glaciale across the Atlantic over distances up to 3940 km (Quintela et al., 2024). Likewise, no genetic structure at mitochondrial cytochrome b was displayed in the North Pacific Ocean in the myctophids Diaphus thetaStenobrachius leucopsarusS. nannochir displaying different diel vertical migration patterns (Kojima et al., 2009). Lack of microsatellite genetic differentiation coupled with large genetic diversity was also reported for the southern hemisphere myctophid Electrona antarctica (Van de Putte et al., 2012). The homogenizing effect of the Southern Coastal Current in the Southern Ocean not only affects this species, but also others such as humped rockcod Gobionotothen gibberifrons (Matschiner et al., 2009) or Antarctic silverfish Pleuragramma antarctica (Zane et al., 2006) for which gene flow mediated by larval dispersal resulted in weak or absent genetic structure. In M. muelleri, the only northern deviating sample was Vesterålen, which seemed to consist of a physical mixture of genetically distinct fjord and offshore individuals. It could be hypothesized that individuals overcoming the fjord sills could be drifted northwards by the Norwegian Current and retained in Vesterålen area, eventually in combination with the effect of the intensive upwelling that brings deep-see water to Lofoten and Vesterålen (Falk & Nøst, 2013). Otherwise, the differentiation detected in the Atlantic Ocean restricts to the Eastern Southern samples and is driven by the influence of the Mediterranean Sea.

      Very little genetic differentiation was observed among samples from the Norwegian fjords, where extremely low but significant differentiation was only detected between the two most distant sites (Osterfjord vs. Boknafjord, FST=0.007**). Our result is in agreement with data from an earlier study by Watkins et al. (1996) who reported a lack of genetic structure in M. muelleri in the western Norwegian fjords using allozymes. Also with allozymes, Kristoffersen and Salvanes (2009) did not detect differentiation for B. glaciale among the fjords, which contrasts with nearly half of the pairwise fjord comparisons revealing significant structure when genotyped with SNPs (Quintela et al., 2024) and with Suneetha and Salvanes (2001) where genetic homogeneity was only detected in fjords with sill depths exceeding 130 m. The sills, typically present at the mouths of the fjords that hinder deep-water exchange with the adjoining coastal areas, could be acting as a physical barrier limiting gene flow within B. glaciale (Kristoffersen & Salvanes, 2009), however, they do not seem to have any effect on M. muelleri. Differential swimming capacity of both species could be partially responsible for this difference in genetic structure. Tracking of individuals in deep water showed that B. glaciale was conspicuously inactive and drifted back and forth with weak tidal currents, essentially acting as plankton while active swimming occasionally occurred in the vertical direction with speeds ranging between <0.5–1 body length s−1 (Kaartvedt et al., 2009). In contrast, M. muelleri swimming capacity has been reported to reach speeds of 2–7 body length·s-1 (i.e. 8–30 cm·s-1) (Torgersen & Kaartvedt, 2001). In addition to larger swimming capacity, M. muelleri has been reported to stay in upper levels in the water column than B. glaciale both during day and night time, not only in the Norwegian fjords but also in the Mediterranean and the Atlantic coast (Giske et al., 1990Staby et al., 2011Bañón et al., 2016Rábade Uberos et al., 2021Olivar et al., 2022Bernal et al., 2023Kapelonis et al., 2023) even reaching the surface in the Norwegian fjords during night time (Staby et al., 2011). Furthermore, M. muelleri has the capability of changing its behaviour in response to ontogeny and internal state (satiation and hunger) as well as to external stimuli (Staby et al., 2011).

      In addition, or maybe because of the reported genetic differentiation across habitats, life history of M. muelleri is known to differ geographically, with maximum age reached in Norwegian waters (Gjøsæter & Kawaguchi, 1980). Along the Norwegian coast, the hydrological and topographic characteristics of the fjords may create unique habitat conditions (Farmer & Freeland, 1983) and differentiation between offshore and fjords was detected in traits such as sexual dimorphism (Kristoffersen & Salvanes, 2001), maturation stage, gonadal investment, relative fecundity and growth rate, which displayed increasing trends from the Norwegian Sea to the coast and towards the fjords (Salvanes & Stockley, 1996). The reproductive investment was higher in a fjord with higher predation risk (Bjelland, 1995) and fjord fish also matured earlier in the season than offshore (Salvanes & Stockley, 1996). Mortality rate was reported higher in offshore areas than within the fjords, which was attributed to fish occupying a depth with higher light levels and therefore more vulnerable to predators (Kristoffersen & Salvanes, 1998). Differences in growth, condition and gonad weight indicated different resource levels caused by different population densities (Kristoffersen & Salvanes, 1998). However, back in the day, life history differences were suggested to be phenotypic and not genetically-driven due to the lack of variation fjord-offshore displayed by allozymes (Watkins et al., 1996), which contrasts with the significant differentiation detected with SNPs in this study (average FST≈0.09). Interestingly, in the same paper, Watkins et al. (1996) detected fjord-offshore differentiation in B. glaciale later confirmed by SNPs (Quintela et al., 2024). Fjord vs. off-shore differentiation has also been described in other marine taxa with high dispersal capacity such as, e.g., tunicate Ciona intestinalis (Johannesson et al., 2018), northern shrimp Pandalus borealis (Hansen et al., 2021), haddock Melanogrammus aeglefinus (Berg et al., 2021), Atlantic cod (Ruzzante et al., 1997Westgaard & Fevolden, 2007Pampoulie et al., 2011) or European sprat (Quintela et al., 2020Pettersson et al., 2024).

      Genetic differentiation across areas also mirrors the differentiation found in life history traits in M. muelleri from the Bay of Biscay, Celtic Sea, and around 59–63 °N in the Norwegian Sea where length-weight relationships revealed differences associated with the fish’s origin, paralleling the annual and daily otolith growth. Likewise, von Bertalanffy growth models parameters increased progressively northwards, in accordance with Bergmann’s rule (Alvarez et al., 2024). The authors suggested the existence of separated units, either genetically or morphologically, representing differences in biological parameters as a signal of geographical divergence. Our genetic data confirms that the Bay of Biscay significantly differed from the Celtic Sea as well as from Norwegian coastal samples within the same latitudinal range; likewise, the Celtic Sea differed from one of the Norwegian samples.

      Oceanographic features like current patterns and discontinuities shape the genetic connectivity in the Mediterranean Sea (Galarza et al., 2009aSchunter et al., 2011) with the Almería-Oran Front (AOF) being the major point of genetic break between the Atlantic Ocean and the Mediterranean (Patarnello et al., 2007). In addition, the entry of less saline water from the Atlantic through the shallow Strait of Gibraltar represents a barrier to gene flow for multiple species (Galarza et al., 2009aGalarza et al., 2009bMarie et al., 2016) and could account for the much lower mesopelagic fish diversity found in the Mediterranean compared with the adjacent Atlantic waters (Olivar et al., 2022). Along the Greek coastline, the Ionian and Aegean Seas shape a complex ecosystem combining a highly irregular coastline and semi-isolated deep basins leading to local differences in species composition (Somarakis et al., 2011Kapelonis et al., 2023). Genetic differentiation within species has thus been attributed to a combination of historic demographic processes as well as hydrological and ecological traits (see Sarropoulou et al., 2022 and references therein). M. muelleri is not oblivious to these factors, which drive the large and highly significant differentiation detected across the Mediterranean (FST ranging from 0.05 to 0.18) except for the pair N. Aegean Sea and Cretan Sea, the only samples that do not seem to experience conspicuous physical barriers to gene flow (see Fig. 1). The Gulf of Corinth (Ionian Sea) significantly differed from all the rest (FST ranging from 0.05 to 0.06) as it is semi-closed body water with limited communication to the open sea via shallow and narrow channels, which configurates unique hydrological and topographical traits (Drakopoulos & Lascaratos, 1998Ramfos et al., 2005). The isolation of this site has formerly revealed in both M. muelleri and B. glaciale by high number of private mtDNA alleles traits (Sarropoulou et al., 2022), in nuclear markers in B. glaciale (Quintela et al., 2024), and in the small-spotted catshark Scyliorhinus canicula (Kousteni et al., 2015). Among all Mediterranean samples, genetic connectivity seems extremely hampered towards the N Euboean Gulf (FST ranging from 0.16 to 0.18), sample that is also genetically closer to the fjord and Atlantic ones.

      Outlier detection and environmental association analyses

      The identification of 20 SNPs as candidate outliers to positive selection, alongside 25 candidates to balancing selection highlights the complexity of selective pressures acting on different loci. Adaptive traits relevant to specific environmental contexts are likely maintained through positive selection whereas loci under balancing selection promote genetic diversity, possibly enabling populations to adapt to fluctuating environmental conditions. This distinction emphasizes the necessity of targeting diverse loci when studying adaptation, as it underscores the multifaceted nature of evolutionary pressures in different habitats.

      The neutral markers retained the patterns of habitat differentiation observed in the total dataset thus highlighting their role of in shaping genetic structure across habitats. Although PCA and DAPC plots revealed differentiation among habitats, their limitations in capturing the Atlantic-Mediterranean transition indicate a pivotal role of historical factors, such as gene flow and connectivity. This suggests that neutral genetic markers can provide valuable but partial insights into the evolutionary history of the populations. The reduced discrimination capacity signals that more nuanced factors beyond simple neutrality may shape genetic differentiation and warrants further investigation into historical influences.

      The candidate loci to positive selection did not have the ability to discriminate between fjords and oceanic samples from the Northern Atlantic cluster. This could be explained by 95% of the candidate outliers being related to temperature according to LFMM and showing extremely similar allele frequencies in the geographic range covering from Vesterålen southwards to the Celtic Sea and westwards to Flemish Cap including both fjord and ocean samples. Likewise, allele frequency clines revealed in 42% of the loci analysed in this study suggest variations due to environmental factors along the latitudinal gradient. In addition to loci linked to temperature, RDA identified a locus seemingly linked to chlorophyll concentration, in agreement with Alamanellis-Zisimopoulos (2023) who reported that bottom depth was the variable with largest importance, followed by the surface chlorophyll concentration (mg/m3) and sea current geostrophic velocity (m/s) when modelling the potential habitat of the species in the Greek seas.

      An interesting feature is the deviating allele frequencies found in the Greek sample from the N. Euboean Gulf, which tended towards the range of values found in the fjords-Atlantic in contrast with the remaining Mediterranean samples. This pattern occurred in some of the candidate loci to positive selection/linked to temperature while being absent in the neutral ones and emphasizes the potential importance of environmental adaptation. The N. Euboean Gulf is a ∼450 m deep basin that connects with adjacent seas through shallow channels and host only two mesopelagic fish species (M. muelleri and B. glaciale). The temperature is some ≈2 °C lower than in other Mediterranean regions and it is the only Greek basin that features hypoxic conditions at depths deeper than 350 m (see Suppl. Fig. S12). The lower temperature in combination with the hypoxic conditions and the geographic isolation of the gulf could account for the observed patterns.

      Geographic distance explained from 58% (neutral loci) to 61% (all markers) of the genetic differentiation under an Isolation-by-Distance (IBD) pattern for the full set of samples. Likewise, temperature correlated with genetic differentiation driving an Isolation-by-Environment (IBE) pattern both at neutral loci, candidate loci to positive selection and all loci (Suppl. Table S6). However, the correlation between geographic distance and temperature (rxy=0.79, P<0.001) makes it difficult to disentangle the effect of both variables. Partial Mantel tests correcting IBD by temperature were significant for all loci, whereas when correcting IBE by geographic distance no correlation whatsoever was found at neutral markers suggesting that the relationship between temperature and genetic differentiation disappears when accounting for geographic distance. Candidate outliers to positive selection retained significant when correcting by geographic distance, what highlights the effect of temperature in habitat-driven differentiation in M. muelleri.

      Genetic response to environmental changes can happen at a high pace in the marine realm. One century was enough to generate a distinct population of European sprat after colonizing an artificially created brackish lake (Quintela et al., 2021Pettersson et al., 2024). The remarkable response of M. muelleri to variable thermal environments suggests that understanding these dynamics may help predict how the species will respond to future climate scenarios and highlights that local adaptation will be critical in informing conservation strategies in the face of ongoing climate change.

      Management implications

      Although strongly debated (Pauly et al., 2021), the ambitious estimates of 10 Gt of mesopelagic fish biomass, which is equivalent to ∼100 times the annual catch of all existing fisheries (Hidalgo & Browman, 2019) fuels the interest of harvesting this relatively intact putative resource, mainly as animal feed for aquaculture (Alvheim et al., 2020Grimaldo et al., 2020). Thus, trial fisheries for mesopelagic fish (including M. muelleri) have recently been conducted by Norwegian vessels in the North Atlantic (Standal & Grimaldo, 2021). In the face of an eventual commercial exploitation, and to ensure sustainability, a compelling body of literature highlights the importance of accurately aligning genetic and biological units when defining fisheries stocks (e.g. Allendorf et al., 2008Hauser & Carvalho, 2008Waples et al., 2008Reiss et al., 2009Funk et al., 2012Cadrin et al., 2014Bernatchez et al., 2017Kerr et al., 2017). The finding of three distinct genetic populations of M. muelleri in the North Atlantic and Mediterranean basins provides background information while suggesting that demographic properties should also be mapped to provide a more comprehensive stock outline. Mesopelagic fishes provide multiple ecosystem services, such as regulation ecosystem services, due to their role in the Biological Carbon Pump and carbon sequestration (John et al., 2016), or supporting ecosystem services, since they are prey for a diversity of predators, including economically important species (Iglesias et al., 2023). In consequence, the ecological importance of M. muelleri (as of many other mesopelagic fish species) must be thoroughly addressed and put into perspective before attempting any harvest with potential harmful consequences.

      Supporting information

      Supplement[supplements/639266_file10.docx]

      DATA AVAILABILITY

      The genotype raw data used in this study can be publicly accessed from the electronic archive of the Institute of Marine Research upon acceptance.

      BENEFIT-SHARING STATEMENT

      Benefits generated: A research collaboration was developed with all the scientists from the countries providing genetic samples, all collaborators are included as co-authors, the results of research have been shared with the provider communities and the broader scientific community, and the research addresses a priority concern, in this case the conservation of organisms being studied. More broadly, our group is committed to international scientific partnerships, as well as institutional capacity building.

      CONFLICTS OF INTEREST STATEMENT

      The authors declare no conflicts of interest.

      ACKNOWLEDGEMENTS

      This study was primarily funded by the Norwegian Department for Trade and Fisheries, and further supported by the EU through funding for the MEESO project through the EU H2020 Research and Innovation Programme, Grant Agreement No 817669. Samples from the Western Mediterranean were collected during a survey performed within the SUMMER EU H2020 Research and Innovation Programme, Grant Agreement No 817806. Samples from the Eastern Mediterranean were obtained under the project MesoBED, funded by the Hellenic Foundation for Research and Innovation and the General Secretariat of Research and Innovation (Greece) (Project No. 449). The sample from Moroccan waters was collected through the scientific surveys with the research vessel Dr. Fridtjof Nansen as part of the collaboration between the Food and Agriculture Organization of the United Nations (FAO) on behalf of the EAF-Nansen Programme and the Kingdom of Morocco. The EAF-Nansen Programme is a partnership between the FAO, the Norwegian Agency for Development Cooperation (Norad), and the Institute of Marine Research (IMR) in Norway for sustainable management of the fisheries in partner countries and regions. The sample from the Flemish Cap was collected by CISC (Spain) on behalf of the European Union sampling programme during a survey conducted onboard the research vessel Vizconde de Eza. We are grateful to Nikolaos Nikolioudakis (IMR) and Iñaki Mendibil (AZTI) for their valuable help with sampling logistics.

      REFERENCES

      Zane, L., Marcato, S., Bargelloni, L., Bortolotto, E., Papetti, C., Simonato, M., … Patarnello, T. (2006). Demographic history and population structure of the Antarctic silverfish Pleuragramma antarcticum. Molecular Ecology, 15(14), 4499–4511. doi:10.1111/j.1365-294X.2006.03105.x

      Alamanellis-Zisimopoulos, A. (2023). Modelling the distribution of the mesopelagic fish Maurolicus muelleri in the Greek seas. (Msc Msc “Oceanography and Management of the Marine Environment”), National and Kapodistrian University of Athens, Athens.

      Allendorf, F. W., England, P. R., Luikart, G., Ritchie, P. A., & Ryman, N. (2008). Genetic effects of harvest on wild animal populations. Trends in Ecology & Evolution, 23(6), 327–337. doi:10.1016/j.tree.2008.02.008

      Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Basic local alignment search tool. Journal of Molecular Biology, 215(3), 403–410. doi:10.1016/S0022-2836(05)80360-2

      Alvarez, P., Aldanondo, N., Wieczorek, A. M., Cariou, T., Boyra, G., Grimaldo, E., … Klevjer, T. (2024). Exploring growth patterns of Maurolicus muelleri across three Northeast Atlantic regions. Fishes, 9(7), 250. doi:10.3390/fishes9070250

      Alvheim, A. R., Kjellevold, M., Strand, E., Sanden, M., & Wiech, M. (2020). Mesopelagic species and their potential contribution to food and feed security—A case study from Norway. Foods, 9(3). Retrieved from doi:10.3390/foods9030344

      Assis, J., Tyberghein, L., Bosch, S., Verbruggen, H., Serrão, E. A., & De Clerck, O. (2018). Bio-ORACLE v2.0: Extending marine data layers for bioclimatic modelling. Global Ecology and Biogeography, 27(3), 277–284. doi:10.1111/geb.12693

      Baliño, B. (1993). Winter distribution and migration of the sound scattering layers, zooplankton and micronekton in Masfjorden, western Norway. Marine Ecology Progress Series, 1+2, 35-50.

      Bañón, R., Arronte, J. C., Rodríguez-Cabello, C., Piñeiro, C. G., Punzón, A., & Serrano, A. (2016). Commented checklist of marine fishes from the Galicia Bank seamount (NW Spain). Zootaxa, 4067(3), 293–333. doi: 10.11646/zootaxa.4067.3.2

      Bargelloni, L., Alarcon, J. A., Alvarez, M. C., Penzo, E., Magoulas, A., Reis, C., & Patarnello, T. (2003). Discord in the family Sparidae (Teleostei): divergent phylogeographical patterns across the Atlantic– Mediterranean divide. Journal of Evolutionary Biology, 16(6), 1149–1158. doi:10.1046/j.1420-9101.2003.00620.x

      Barros-García, D., Froufe, E., Bañón, R., Arronte, J. C., Baldó, F., & de Carlos, A. (2020). Phylogeography highlights two different Atlantic/Mediterranean lineages and a phenotypic latitudinal gradient for the deep-sea morid codling Lepidion lepidion (Gadiformes: Moridae). Deep Sea Research Part I: Oceanographic Research Papers, 157, 103212. doi:10.1016/j.dsr.2019.103212

      Benestan, L. M., Rougemont, Q., Senay, C., Normandeau, E., Parent, E., Rideout, R., Parent, G. J. (2021). Population genomics and history of speciation reveal fishery management gaps in two related redfish species (Sebastes mentella and Sebastes fasciatus). Evolutionary Applications, 14(2), 588–606. doi:10.1111/eva.13143

      Benjamini, Y., & Hochberg, Y. (1995). Controlling the False Discovery Rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Methodological), 57(1), 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x

      Berg, P. R., Jorde, P. E., Glover, K. A., Dahle, G., Taggart, J. B., Korsbrekke, K., Westgaard, J.-I. (2021). Genetic structuring in Atlantic haddock contrasts with current management regimes. ICES Journal of Marine Science, 78(1), 1–13. doi:10.1093/icesjms/fsaa204

      Bernal, A., Tuset, V. M., & Olivar, M. P. (2023). Multiple approaches to the trophic role of mesopelagic fish around the Iberian Peninsula. Animals, 13(5), 886. doi:10.3390/ani13050886

      Bernatchez, L., Wellenreuther, M., Araneda, C., Ashton, D. T., Barth, J. M. I., Beacham, T. D., … Withler, R. E. (2017). Harnessing the power of genomics to secure the future of seafood. Trends in Ecology & Evolution, 32(9), 665–680. doi:10.1016/j.tree.2017.06.010

      Bjelland, O. (1995). Life-history tactics of two fjordic populations of Maurolicus muelleri. (Ms. Master thesis), University of Bergen, Bergen.

      Boyd, P. W. (2015). Toward quantifying the response of the oceans’ biological pump to climate change. Frontiers in Marine Science, 2. doi:10.3389/fmars.2015.00077

      Boyd, P. W., Claustre, H., Levy, M., Siegel, D. A., & Weber, T. (2019). Multi-faceted particle pumps drive carbon sequestration in the ocean. Nature, 568(7752), 327–335. doi:10.1038/s41586-019-1098-2

      Cadrin, S. X., Kerr, L. A., & Mariani, S. (2014). Stock Identification Methods: Applications in Fishery Science (Second edition ed.): Academic Press.

      Chen, S., Zhou, Y., Chen, Y., & Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17), i884–i890. doi:10.1093/bioinformatics/bty560

      Cherel, Y., Fontaine, C., Richard, P., & Labatc, J.-P. (2010). Isotopic niches and trophic levels of myctophid fishes and their predators in the Southern Ocean. Limnology and Oceanography, 55(1), 324–332. doi:10.4319/lo.2010.55.1.0324

      Dahle, G., Quintela, M., Johansen, T., Westgaard, J.-I., Besnier, F., Aglen, A., … Glover, K. A. (2018). Analysis of coastal cod (Gadus morhua L.) sampled on spawning sites reveals a genetic gradient throughout Norway’s coastline. BMC Genetics, 19(1), 42. doi:10.1186/s12863-018-0625-8

      Drakopoulos, P. G., & Lascaratos, A. (1998). A preliminary study on the internal tides of the gulfs of Patras and Korinthos, Greece. Continental Shelf Research, 18(12), 1517–1529. doi:10.1016/S0278-4343(98)00052-1

      Dray, S., & Dufour, A.-B. (2007). The ade4 Package: Implementing the Duality Diagram for Ecologists. Journal of Statistical Software, 22(4), 1–20. doi:10.18637/jss.v022.i04

      Drazen, J. C., & Sutton, T. T. (2017). Dining in the deep: The feeding ecology of deep-sea fishes. Annual Review of Marine Science, 9, 337–366. doi:10.1146/annurev-marine-010816-060543

      Ducklow, H. (2001). Upper ocean carbon export and the Biological Pump. Oceanography, 14, 50–54. doi:10.5670/oceanog.2001.06

      Excoffier, L., Laval, G., & Schneider, S. (2005). Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online, 1, 47–50. doi:10.1177/117693430500100003

      Falk, A. H., & Nøst, O.-A. (2013). Oppstrømming av dyphavsvann – litteraturstudie av oppstrømming utenfor Salten/Lofoten/Vesterålen. Akvaplan-niva rapport nr 6311–01, 32 pp.

      FAO. (2020). Report of the Working Group on the Assessment of Small Pelagic Fish of Northwest Africa Casablanca, Morocco, 8–13 July 2019. Rapport de groupe de travail sur l’évaluation des petits pêlagiques au large de l’Afrique Nord-Occidentale Casablanca, Maroc, 8-13 juillet 2019. Fishery Committee for the Eastern Central Atlantic (CECAF)/Comité des pêches pour l’Atlantique Centre-Est (COPACE). RFAO Fisheries and Aquaculture Report No. 1309/FAO, Rapport sur les pêches et l’aquaculture no 1309. . Retrieved from Rome:

      Farmer, D. M., & Freeland, H. J. (1983). The physical oceanography of Fjords. Progress in Oceanography, 12(2), 147–219. doi:10.1016/0079-6611(83)90004-6

      Flynn, A. J., & Marshall, N. J. (2013). Lanternfish (Myctophidae) zoogeography off Eastern Australia: A comparison with physicochemical biogeography. PLOS ONE, 8(12), e80950. doi:10.1371/journal.pone.0080950

      Foll, M., & Gaggiotti, O. (2006). Identifying the environmental factors that determine the genetic structure of populations. Genetics, 174(2), 875–891. doi:10.1534/genetics.106.059451

      Foll, M., & Gaggiotti, O. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics, 180(2), 977–993. doi:10.1534/genetics.108.092221

      Fontaine, M. C. (2016). Chapter Eleven – Harbour porpoises, Phocoena phocoena, in the Mediterranean Sea and adjacent regions: Biogeographic relicts of the Last Glacial Period. In G. Notarbartolo Di Sciara, M. Podestà, & B. E. Curry (Eds.), Advances in Marine Biology (Vol. 75, pp. 333–358): Academic Press.

      Forester, B. R., Lasky, J. R., Wagner, H. H., & Urban, D. L. (2018). Comparing methods for detecting multilocus adaptation with multivariate genotype–environment associations. Molecular Ecology, 27(9), 2215–2233. doi:10.1111/mec.14584

      Frichot, E., & François, O. (2015). LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution, 6(8), 925–929. doi:10.1111/2041-210X.12382

      Frichot, E., Schoville, S. D., Bouchard, G., & François, O. (2013). Testing for associations between loci and environmental gradients using latent factor mixed models. Molecular biology and evolution, 30(7), 1687–1699. doi:10.1093/molbev/mst063

      Funk, W. C., McKay, J. K., Hohenlohe, P. A., & Allendorf, F. W. (2012). Harnessing genomics for delineating conservation units. Trends in Ecology & Evolution, 27(9), 489–496. doi:10.1016/j.tree.2012.05.012

      Gabriel, S., Ziaugra, L., & Tabbaa, D. (2009). SNP Genotyping Using the Sequenom MassARRAY iPLEX Platform. Current Protocols in Human Genetics, 60(1), 2.12.11–12.12.18. doi:10.1002/0471142905.hg0212s60

      Galarza, J. A., Carreras-Carbonell, J., Macpherson, E., Pascual, M., Roques, S., Turner, G. F., & Rico, C. (2009a). The influence of oceanographic fronts and early-life-history traits on connectivity among littoral fish species. Proceedings of the National Academy of Sciences, 106(5), 1473–1478. doi:10.1073/pnas.0806804106

      Galarza, J. A., Turner, G. F., Macpherson, E., & Rico, C. (2009b). Patterns of genetic differentiation between two co-occurring demersal species: the red mullet (Mullus barbatus) and the striped red mullet (Mullus surmuletus). Canadian Journal of Fisheries and Aquatic Sciences, 66(9), 1478–1490. doi:10.1139/F09-098

      Giske, J., Aksnes, D. L., Baliño, B. M., Kaartvedt, S., Lie, U., Nordeide, J. T., … Aadnesen, A. (1990). Vertical distribution and trophic interactions of zooplankton and fish in Masfjorden, Norway. Sarsia, 75, 65–81.

      Gjøsæter, J., & Kawaguchi, K. (1980). A review of the world resources of mesopelagic fish. FAO Fish. Tech. Pap. No. 193. FIRM/TI93. 151 pp.

      Goodson, M. S., Giske, J., & Rosland, R. (1995). Growth and ovarian development of Maurolicus muelleri during spring. Marine Biology, 124(2), 185–195. doi: 10.1007/BF00347122

      Grimaldo, E., Grimsmo, L., Álvarez, P., Herrmann, B., Møen Tveit, G., Tiller, R., … Selnes, M. (2020). Investigating the potential for a commercial fishery in the Northeast Atlantic utilizing mesopelagic species. ICES Journal of Marine Science, 77(7-8), 2541–2556. doi:10.1093/icesjms/fsaa114

      Hansen, A., Westgaard, J.-I., Søvik, G., Hanebrekke, T., Nilssen, E. M., Jorde, P. E., … Johansen, T. (2021). Genetic differentiation between inshore and offshore populations of northern shrimp (Pandalus borealis). ICES Journal of Marine Science, 78(9), 3135–3146. doi:10.1093/icesjms/fsab181

      Hauser, L., & Carvalho, G. R. (2008). Paradigm shifts in marine fisheries genetics: ugly hypotheses slain by beautiful facts. Fish and Fisheries, 9(4), 333–362. doi:10.1111/j.1467-2979.2008.00299.x

      Hays, G. C. (2003). A review of the adaptive significance and ecosystem consequences of zooplankton diel vertical migrations. Hydrobiologia, 503(1), 163–170. doi:10.1023/B:HYDR.0000008476.23617.b0

      Hidaka, K., Kawaguchi, K., Murakami, M., & Takahashi, M. (2001). Downward transport of organic carbon by diel migratory micronekton in the western equatorial Pacific: its quantitative and qualitative importance. Deep Sea Research Part I: Oceanographic Research Papers, 48(8), 1923–1939. doi:10.1016/S0967-0637(01)00003-6

      Hidalgo, M., & Browman, H. I. (2019). Developing the knowledge base needed to sustainably manage mesopelagic resources. ICES Journal of Marine Science, 76(3), 609–615. doi:10.1093/icesjms/fsz067

      Iglesias, I. S., Santora, J. A., Fiechter, J., & Field, J. C. (2023). Mesopelagic fishes are important prey for a diversity of predators. Frontiers in Marine Science, 10. doi:10.3389/fmars.2023.1220088

      Irigoien, X., Klevjer, T. A., Røstad, A., Martinez, U., Boyra, G., Acuña, J. L., … Kaartvedt, S. (2014). Large mesopelagic fishes biomass and trophic efficiency in the open ocean. Nature Communications, 5(1), 3271. doi:10.1038/ncomms4271

      Jakobsson, M., & Rosenberg, N. A. (2007). CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23(14), 1801– 1806. doi:10.1093/bioinformatics/btm233

      Jenkins, T. L., Ellis, C. D., Triantafyllidis, A., & Stevens, J. R. (2019). Single nucleotide polymorphisms reveal a genetic cline across the north-east Atlantic and enable powerful population assignment in the European lobster. Evolutionary Applications, 12(10), 1881–1899. doi:10.1111/eva.12849

      Johannesson, K., & Andre, C. (2006). Life on the margin: genetic isolation and diversity loss in a peripheral marine ecosystem, the Baltic Sea. Molecular Ecology, 15(8), 2013–2029.

      Johannesson, K., Ring, A.-K., Johannesson, K. B., Renborg, E., Jonsson, P. R., & Havenhand, J. N. (2018). Oceanographic barriers to gene flow promote genetic subdivision of the tunicate Ciona intestinalis in a North Sea archipelago. Marine Biology, 165(8), 126. doi:10.1007/s00227-018-3388-x

      John, M., Borja, Á., Chust, G., Heath, M., Grigorov, I., Mariani, P., … Santos, R. (2016). A dark hole in our understanding of marine ecosystems and their services: Perspectives from the mesopelagic community. Frontiers in Marine Science, 3. doi:10.3389/fmars.2016.00031

      Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403–1405. doi:10.1093/bioinformatics/btn129

      Jombart, T., & Collins, C. (2015). A tutorial for discriminant analysis of principal components (DAPC) using adegenet 2.0.0https://adegenet.r-forge.r-project.org/files/tutorial-dapc.pdf.

      Jombart, T., Devillard, S., & Balloux, F. (2010). Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics, 11(1), 94. doi:10.1186/1471-2156-11-94

      Jorde, P. E., Huserbråten, M. B. O., Seliussen, B. B., Myksvoll, M. S., Vikebø, F. B., Dahle, G., … Johansen, T. (2021). The making of a genetic cline: introgression of oceanic genes into coastal cod populations in the Northeast Atlantic. Canadian Journal of Fisheries and Aquatic Sciences, 78(7), 958–968. doi:10.1139/cjfas-2020-0380

      Kamvar, Z. N., Tabima, J. F., & Grünwald, N. J. (2014). Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ, 2, e281. doi:10.7717/peerj.281

      Kapelonis, Z., Siapatis, A., Machias, A., Somarakis, S., Markakis, K., Giannoulaki, M., … Tsagarakis, K. (2023). Seasonal patterns in the mesopelagic fish community and associated deep scattering layers of an enclosed deep basin. Scientific Reports, 13(1), 17890. doi:10.1038/s41598-023-44765-5

      Kerr, L. A., Hintzen, N. T., Cadrin, S. X., Clausen, L. W., Dickey-Collas, M., Goethel, D. R., … Nash, R. D. M. (2017). Lessons learned from practical approaches to reconcile mismatches between biological population structure and stock units of marine fish. ICES Journal of Marine Science, 74(6), 1708–1722. doi:10.1093/icesjms/fsw188

      Knaus, B. J., & Grünwald, N. J. (2016). VcfR: a package to manipulate and visualize VCF format data in R. bioRxiv, 041277. doi:10.1101/041277

      Knaus, B. J., & Grünwald, N. J. (2017). vcfr: a package to manipulate and visualize variant call format data in R. Mol Ecol Resour, 17(1), 44–53. doi:10.1111/1755-0998.12549

      Knutsen, H., Jorde, P. E., Hutchings, J. A., Hemmer-Hansen, J., Grønkjær, P., Jørgensen, K.-E. M., … Olsen, E. M. (2018). Stable coexistence of genetically divergent Atlantic cod ecotypes at multiple spatial scales. Evolutionary Applications, 11(9), 1527–1539. doi:10.1111/eva.12640

      Knutsen, H., Jorde, P. E., Sannæs, H., Rus Hoelzel, A., Bergstad, O. A., Stefanni, S., … Stenseth, N. C. (2009). Bathymetric barriers promoting genetic structure in the deepwater demersal fish tusk (Brosme brosme). Molecular Ecology, 18(15), 3151–3162. doi:10.1111/j.1365-294X.2009.04253.x

      Kojima, S., Moku, M., & Kawaguchi, K. (2009). Genetic diversity and population structure of three dominant myctophid fishes (Diaphus theta, Stenobrachius leucopsarus, and S. nannochir) in the North Pacific Ocean. Journal of Oceanography, 65(2), 187–193. doi:10.1007/s10872-009-0018-8

      Koubbi, P., Moteki, M., Duhamel, G., Goarant, A., Hulley, P.-A., O’Driscoll, R., … Hosie, G. (2011). Ecoregionalization of myctophid fish in the Indian sector of the Southern Ocean: Results from generalized dissimilarity models. Deep Sea Research Part II: Topical Studies in Oceanography, 58(1), 170–180. doi:10.1016/j.dsr2.2010.09.007

      Kousteni, V., Kasapidis, P., Kotoulas, G., & Megalofonou, P. (2015). Strong population genetic structure and contrasting demographic histories for the small-spotted catshark (Scyliorhinus canicula) in the Mediterranean Sea. Heredity, 114(3), 333–343. doi:10.1038/hdy.2014.107

      Kristoffersen, J. B., & Salvanes, A. G. V. (1998). Life history of Maurolicus muelleri in fjordic and oceanic environments. Journal of Fish Biology, 53(6), 1324–1341. doi:10.1111/j.1095-8649.1998.tb00252.x

      Kristoffersen, J. B., & Salvanes, A. G. V. (2001). Sexual size dimorphism and sex ratio in Müller’s pearlside (Maurolicus muelleri). Marine Biology, 138(6), 1087–1092. doi:10.1007/s002270000529

      Kristoffersen, J. B., & Salvanes, A. G. V. (2009). Distribution, growth, and population genetics of the glacier lanternfish (Benthosema glaciale) in Norwegian waters: Contrasting patterns in fjords and the ocean. Marine Biology Research, 5(6), 596–604. doi:10.1080/17451000903042479

      Kaartvedt, S., Røstad, A., Klevjer, T., & Staby, A. (2009). Use of bottom-mounted echo sounders in exploring behavior of mesopelagic fishes. Marine Ecology Progress Series, 395, 109–118. doi:10.3354/meps08174

      Lapègue, S., Reisser, C., Harrang, E., Heurtebise, S., & Bierne, N. (2023). Genetic parallelism between European flat oyster populations at the edge of their natural range. Evolutionary Applications, 16(2), 393–407. doi:10.1111/eva.13449

      Le Moan, A., Gagnaire, P.-A., & Bonhomme, F. (2016). Parallel genetic divergence among coastal–marine ecotype pairs of European anchovy explained by differential introgression after secondary contact. Molecular Ecology, 25(13), 3187–3202. doi:10.1111/mec.13627

      Li, D., Liu, C.-M., Luo, R., Sadakane, K., & Lam, T.-W. (2015). MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics, 31(10), 1674–1676. doi:10.1093/bioinformatics/btv033

      Li, D., Luo, R., Liu, C. M., Leung, C. M., Ting, H. F., Sadakane, K., … Lam, T. W. (2016). MEGAHIT v1.0: A fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods, 102, 3–11. doi:10.1016/j.ymeth.2016.02.020

      Li, H. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics, 27(21), 2987–2993. doi:10.1093/bioinformatics/btr509

      Li, H., & Durbin, R. (2010). Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics, 26(5), 589–595. doi:10.1093/bioinformatics/btp698

      Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., … Genome Project Data Processing, S. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078-2079. doi:10.1093/bioinformatics/btp352

      Li, Y.-L., & Liu, J.-X. (2018). StructureSelector: A web-based software to select and visualize the optimal number of clusters using multiple methods. Molecular Ecology Resources, 18(1), 176–177. doi:10.1111/1755-0998.12719

      Lischer, H. E. L., & Excoffier, L. (2012). PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics, 28(2), 298–299. doi:10.1093/bioinformatics/btr642

      Mantel, N. (1967). The detection of disease of clustering and a generalized regression approach. Cancer Research, 27(2), 209–220.

      Marie, A. D., Lejeusne, C., Karapatsiou, E., Cuesta, J. A., Drake, P., Macpherson, E., … Rico, C. (2016). Implications for management and conservation of the population genetic structure of the wedge clam Donax trunculus across two biogeographic boundaries. Scientific Reports, 6(1), 39152. doi:10.1038/srep39152

      Martin, A., Boyd, P., Buesseler, K., Cetinic, I., Claustre, H., Giering, S., … Guidi, L. (2020). The oceans’ twilight zone must be studied now, before it is too late. Nature, 580(7801), 26–28. doi:10.1038/d41586-020-00915-7

      Matschiner, M., Hanel, R., & Salzburger, W. (2009). Gene flow by larval dispersal in the Antarctic notothenioid fish Gobionotothen gibberifrons. Molecular Ecology, 18(12), 2574–2587. doi:10.1111/j.1365-294X.2009.04220.x

      Meyer, L., Barry, P., Riquet, F., Foote, A., Der Sarkissian, C., Cunha, R. L., Gagnaire, P. A. (2024). Divergence and gene flow history at two large chromosomal inversions underlying ecotype differentiation in the long-snouted seahorse. Mol Ecol, e17277. doi:10.1111/mec.17277

      Miller, J. M., Cullingham, C. I., & Peery, R. M. (2020). The influence of a priori grouping on inference of genetic clusters: simulation study and literature review of the DAPC method. Heredity, 125(5), 269–280. doi:10.1038/s41437-020-0348-2

      Nowicki, M., DeVries, T., & Siegel, D. A. (2022). Quantifying the Carbon export and sequestration pathways of the Ocean’s Biological Carbon Pump. Global Biogeochemical Cycles, 36(3), e2021GB007083. doi:10.1029/2021GB007083

      Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Wagner, H. (2020). vegan community ecology package version 2.5-7. https://CRAN.R-project.org/package=vegan.

      Olivar, M. P., Castellón, A., Sabatés, A., Sarmiento-Lezcano, A., Emelianov, M., Bernal, A., Brierley, A. S. (2022). Variation in mesopelagic fish community composition and structure between Mediterranean and Atlantic waters around the Iberian Peninsula. Frontiers in Marine Science, 9. doi:10.3389/fmars.2022.1028717

      Pampoulie, C., Daníelsdóttir, A. K., Storr-Paulsen, M., Hovgård, H., Hjörleifsson, E., & Steinarsson, B. Æ. (2011). Neutral and nonneutral genetic markers revealed the presence of inshore and offshore stock components of Atlantic cod in Greenland waters. Transactions of the American Fisheries Society, 140(2), 307–319. doi:10.1080/00028487.2011.567850

      Pante, E., & Simon-Bouhet, B. (2013). marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. PloS ONE, 8(9), e73051. doi: 10.1371/journal.pone.0073051

      Patarnello, T., Volckaert, F. A. M. J., & Castilho, R. (2007). Pillars of Hercules: is the Atlantic–Mediterranean transition a phylogeographical break? Molecular Ecology, 16(21), 4426–4444. doi:10.1111/j.1365-294X.2007.03477.x

      Pauly, D., Piroddi, C., Hood, L., Bailly, N., Chu, E., Lam, V., … Palomares, M. L. D. (2021). The biology of mesopelagic fishes and their catches (1950–2018) by commercial and experimental fisheries. Journal of Marine Science and Engineering, 9(10). doi:10.3390/jmse9101057

      Peakall, R., & Smouse, P. E. (2006). GenAlEx 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6(1), 288–295. doi: 10.1111/j.1471-8286.2005.01155.x

      Pereira, P., Teixeira, J., & Velo-Antón, G. (2018). Allele surfing shaped the genetic structure of the European pond turtle via colonization and population expansion across the Iberian Peninsula from Africa. Journal of Biogeography, 45(9), 2202–2215. doi:10.1111/jbi.13412

      Pérez-Losada, M., Guerra, A., Carvalho, G. R., Sanjuan, A., & Shaw, P. W. (2002). Extensive population subdivision of the cuttlefish Sepia officinalis (Mollusca: Cephalopoda) around the Iberian Peninsula indicated by microsatellite DNA variation. Heredity, 89(6), 417–424. doi:10.1038/sj.hdy.6800160

      Pettersson, M. E., Quintela, M., Besnier, F., Deng, Q., Berg, F., Kvamme, C., Andersson, L. (2024). Limited parallelism in genetic adaptation to brackish water bodies in European sprat and Atlantic herring. Genome Biology and Evolution 16(7), evae133. doi:10.1093/gbe/evae133

      Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics, 155(2), 945–959. doi:10.1093/genetics/155.2.945

      Puechmaille, S. J. (2016). The program structure does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. Molecular Ecology Resources, 16(3), 608–627. doi:10.1111/1755-0998.12512

      Quintela, M., García-Seoane, E., Dahle, G., Klevjer, T. A., Melle, W., Lille-Langøy, R., … Glover, K. (2024). Genetics in the ocean’s twilight zone: population structure of the glacier lanternfish Benthosema glaciale Reinhardt, 1837 across its distribution range. Evolutionary Applications.

      Quintela, M., Kvamme, C., Bekkevold, D., Nash, R. D. M., Jansson, E., Sørvik, A. G., Glover, K. A. (2020). Genetic analysis redraws the management boundaries for the European sprat. Evolutionary Applications, 13, 1906–1922. doi:10.1111/eva.12942

      Quintela, M., Richter-Boix, À., Bekkevold, D., Kvamme, C., Berg, F., Jansson, E., Glover, K. A. (2021). Genetic response to human-induced habitat changes in the marine environment: A century of evolution of European sprat in Landvikvannet, Norway. Ecology and Evolution, 11(4), 1691–1718. doi:10.1002/ece3.7160

      Rábade Uberos, S., Vergara Castaño, A. R., Domínguez-Petit, R., & Saborido-Rey, F. (2021). Larval fish community in the Northwestern Iberian Upwelling System during the summer period. Oceans, 2(4), 700–722. doi:10.3390/oceans2040040

      Ramfos, A., Somarakis, S., Koutsikopoulos, C., & Fragopoulu, N. (2005). Summer mesozooplankton distribution in coastal waters of central Greece (eastern Mediterranean). I. Hydrology and group composition. Journal of the Marine Biological Association of the United Kingdom, 85(4), 755–764. doi:10.1017/S0025315405011665

      Rasmussen, O. I., & Giske, J. (1994). Life-history parameters and vertical distribution of Maurolicus muelleri in Masfjorden in summer. Marine Biology, 120(4), 649–664. doi:10.1007/BF00350086

      Rees, D. J., Byrkjedal, I., & Sutton, T. T. (2017). Pruning the Pearlsides: Reconciling morphology and molecules in mesopelagic fishes (Maurolicus: Sternoptychidae). Deep Sea Research Part II: Topical Studies in Oceanography, 137, 246–257. doi:10.1016/j.dsr2.2016.04.024

      Rees, D. J., Poulsen, J. Y., Sutton, T. T., Costa, P. A. S., & Landaeta, M. F. (2020). Global phylogeography suggests extensive eucosmopolitanism in Mesopelagic Fishes (Maurolicus: Sternoptychidae). Scientific Reports, 10(1), 20544. doi:10.1038/s41598-020-77528-7

      Reiss, H., Hoarau, G., Dickey-Collas, M., & Wolff, W. J. (2009). Genetic population structure of marine fish: mismatch between biological and fisheries management units. Fish and Fisheries, 10(4), 361–395. doi:10.1111/j.1467-2979.2008.00324.x

      Rellstab, C., Gugerli, F., Eckert, A. J., Hancock, A. M., & Holderegger, R. (2015). A practical guide to environmental association analysis in landscape genomics. Molecular Ecology, 24(17), 4348–4370. doi:10.1111/mec.13322

      Reuschel, S., Cuesta, J. A., & Schubart, C. D. (2010). Marine biogeographic boundaries and human introduction along the European coast revealed by phylogeography of the prawn Palaemon elegans. Molecular Phylogenetics and Evolution, 55(3), 765–775. doi:10.1016/j.ympev.2010.03.021

      Riesgo, A., Taboada, S., Pérez-Portela, R., Melis, P., Xavier, J. R., Blasco, G., & López-Legentil, S. (2019). Genetic diversity, connectivity and gene flow along the distribution of the emblematic Atlanto-Mediterranean sponge Petrosia ficiformis (Haplosclerida, Demospongiae). BMC Evolutionary Biology, 19(1), 24. doi:10.1186/s12862-018-1343-6

      Riquet, F., Liautard-Haag, C., Woodall, L., Bouza, C., Louisy, P., Hamer, B., … Bierne, N. (2019). Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages. Evolution, 73(4), 817–835. doi:10.1111/evo.13696

      Robinson, C., Steinberg, D. K., Anderson, T. R., Arístegui, J., Carlson, C. A., Frost, J. R., … Zhang, J. (2010). Mesopelagic zone ecology and biogeochemistry – a synthesis. Deep Sea Research Part II: Topical Studies in Oceanography, 57(16), 1504–1518. doi:10.1016/j.dsr2.2010.02.018

      Rodriguez-Ezpeleta, N., Álvarez, P., & Irigoien, X. (2017). Genetic diversity and connectivity in Maurolicus muelleri in the Bay of Biscay inferred from thousands of SNP markers. Frontiers in Genetics, 8. doi:10.3389/fgene.2017.00195

      Rosenberg, M. S., & Anderson, C. D. (2011). PASSaGE: Pattern Analysis, Spatial Statistics and Geographic Exegesis. Version 2. Methods in Ecology and Evolution, 2(3), 229–232. doi: 10.1111/j.2041-210X.2010.00081.x

      Rousset, F. (1997). Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics, 145(4), 1219–1228.

      Rousset, F. (2008). GENEPOP’007: A complete re-implementation of the genepop software for Windows and Linux. Molecular Ecology Resources, 8(1), 103–106. doi: 10.1111/j.1471-8286.2007.01931.x

      Ruzzante, D. E., Taggart, C. T., Cook, D., & Goddard, S. V. (1997). Genetic differentiation between inshore and offshore Atlantic cod (Gadus morhua) off Newfoundland: a test and evidence of temporal stability. Canadian Journal of Fisheries and Aquatic Sciences, 54(11), 2700–2708. doi:10.1139/f97-170

      Salvanes, A. G. V., & Stockley, B. M. (1996). Spatial variation of growth and gonadal developments of Maurolicus muelleri in the Norwegian Sea and in a Norwegian fjord. Marine Biology, 126(2), 321–332. doi:10.1007/BF00347456

      Sanders, R., Henson, S. A., Koski, M., De La Rocha, C. L., Painter, S. C., Poulton, A. J., … Martin, A. P. (2014). The Biological Carbon Pump in the North Atlantic. Progress in Oceanography, 129, 200–218. doi:10.1016/j.pocean.2014.05.005

      Sarropoulou, X., Tsaparis, D., Tsagarakis, K., Badouvas, N., & Tsigenopoulos, C. S. (2022). Different patterns of population structure and genetic diversity of three mesopelagic fishes in the Greek Seas. Mediterranean Marine Science, 23(3), 536–545. doi:10.12681/mms.28567

      Schunter, C., Carreras-Carbonell, J., Macpherson, E., Tintoré, J., Vidal-Vijande, E., Pascual, A., … Pascual, M. (2011). Matching genetics with oceanography: directional gene flow in a Mediterranean fish species. Molecular Ecology, 20(24), 5167–5181. doi:10.1111/j.1365-294X.2011.05355.x

      Shreeve, R. S., Collins, M. A., Tarling, G. A., Main, C. E., Ward, P., & Johnston, N. M. (2009). Feeding ecology of myctophid fishes in the northern Scotia Sea. Marine Ecology Progress Series, 386, 221–236.

      Sigman, D. M., & Haug, G. H. (2006). The Biological Pump in the Past. In H. Elderfield (Ed.), Treatise on Geochemistry. The Oceans and Marine Geochemistry (Vol. 6, pp. 491–528). Cambridge: Elsevier.

      Slatkin, M. (1993). Isolation by Distance in equilibrium and non-equilibrium populations. Evolution, 47(1), 264–279. doi: 10.1111/j.1558-5646.1993.tb01215.x

      Somarakis, S., Isari, S., & Machias, A. (2011). Larval fish assemblages in coastal waters of central Greece: reflections of topographic and oceanographic heterogeneity. Scientia Marina, 75(3), 605–618. doi:10.3989/scimar.2011.75n3605

      Staby, A., Røstad, A., & Kaartvedt, S. (2011). Long-term acoustical observations of the mesopelagic fish Maurolicus muelleri reveal novel and varied vertical migration patterns. Marine Ecology Progress Series, 441, 241–255. doi:10.3354/meps09363

      Standal, D., & Grimaldo, E. (2020). Institutional nuts and bolts for a mesopelagic fishery in Norway. Marine Policy, 119, 104043. doi:10.1016/j.marpol.2020.104043

      Standal, D., & Grimaldo, E. (2021). Lost in translation? Practical- and scientific input to the mesopelagic fisheries discourse. Marine Policy, 134, 104785. doi:10.1016/j.marpol.2021.104785

      Stukel, M. R., Irving, J. P., Kelly, T. B., Ohman, M. D., Fender, C. K., & Yingling, N. (2023). Carbon sequestration by multiple biological pump pathways in a coastal upwelling biome. Nature Communications, 14(1), 2024. doi:10.1038/s41467-023-37771-8

      Suneetha, K. B., & Nævdal, G. (2001). Genetic and morphological stock structure of the pearlside, Maurolicus muelleri (Pisces, Sternoptychidae), among Norwegian fjords and offshore area. Sarsia, 86(3), 191–201. doi:10.1080/00364827.2001.10420475

      Suneetha, K. B., & Salvanes, A. G. V. (2001). Population genetic structure of the glacier lanternfish, Benthosema glaciale (Myctophidae) in Norwegian waters. Sarsia, 86, 203–212.

      Team, R. C. (2020). R: A language and environment for statistical computing. (Version). Vienna: R Foundation for Statistical Computing. Retrieved from https://www.R-project.org/

      Tine, M., Kuhl, H., Gagnaire, P.-A., Louro, B., Desmarais, E., Martins, R. S. T., … Reinhardt, R. (2014). European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation. Nature Communications, 5, 5770. doi:10.1038/ncomms6770 https://www.nature.com/articles/ncomms6770#supplementary-information

      Torgersen, T., & Kaartvedt, S. (2001). In situ swimming behaviour of individual mesopelagic fish studied by split-beam echo target tracking. ICES Journal of Marine Science, 58(1), 346–354. doi:10.1006/jmsc.2000.1016

      Tyberghein, L., Verbruggen, H., Pauly, K., Troupin, C., Mineur, F., & De Clerck, O. (2012). Bio-ORACLE: a global environmental dataset for marine species distribution modelling. Global Ecology and Biogeography, 21(2), 272–281. doi:10.1111/j.1466-8238.2011.00656.x

      Van de Putte, A. P., Van Houdt, J. K. J., Maes, G. E., Hellemans, B., Collins, M. A., & Volckaert, F. A. M. (2012). High genetic diversity and connectivity in a common mesopelagic fish of the Southern Ocean: The myctophid Electrona antarctica. Deep Sea Research Part II: Topical Studies in Oceanography, 59-60, 199–207. doi:10.1016/j.dsr2.2011.05.011

      Waples, R. S., Punt, A. E., & Cope, J. M. (2008). Integrating genetic data into management of marine resources: how can we do it better? Fish and Fisheries, 9(4), 423–449. doi:10.1111/j.1467-2979.2008.00303.x

      Watkins, A., Torkhildsen, S., & Nævdal, G. (1996). Allozyme variation within three species of mesopelagic fishes. Institutt for Fiskeri-og Marinbiologi. Rapport 4, 1–22.

      Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution, 38(6), 1358–1370.

      Westgaard, J.-I., & Fevolden, S.-E. (2007). Atlantic cod (Gadus morhua L.) in inner and outer coastal zones of northern Norway display divergent genetic signature at non-neutral loci. Fisheries Research, 85, 306–315.

      Wiech, M., Silva, M., Meier, S., Tibon, J., Berntssen, M. H. G., Duinker, A., & Sanden, M. (2020). Undesirables in mesopelagic species and implications for food and feed safety-insights from Norwegian fjords. Foods, 9(9). doi:10.3390/foods9091162

      Wright, S. (1943). Isolation by distance. Genetics, 28(2), 114–138. doi:10.1093/genetics/28.2.114