Exploring reported population differences in Norway lobster (Nephrops norvegicus) in the Pomo Pits region of the Adriatic Sea using genome-wide markers

Tissue sampling and DNA extraction

Samples of adult Nephrops were collected from ten sites, including seven from the Adriatic Sea, one from the Aegean Sea, one from the Tyrrhenian Sea, and one from the Firth of Clyde in the northeast Atlantic (Table 1, Figs. 1A–1B). Samples from the Pomo Pits area were collected in November 2016 and October 2019 during the “ScamPo” experimental trawl survey, carried out yearly by CNR IRBIM of Ancona west of the Adriatic midline to monitor the effects of the management measures implemented since 2015. Samples from Ancona and Chioggia were collected in March and July 2019, respectively, by CNR IRBIM staff during biological sampling of commercial catches carried out within the European Data Collection Framework. Tissue samples were obtained by excising two pleopods or one pereiopod. All samples were placed in 95–100% ethanol and stored at 4 °C until DNA extraction. Genomic DNA was extracted using a salting-out protocol (Jenkins, Ellis & Stevens, 2018) and the quality of each extract was assessed on a 0.8% agarose gel. DNA purity was measured using a Nanodrop One and a Qubit 3.0 was used to quantify DNA concentration. The carapace length and sex of each animal sampled was also recorded (except for the Aegean and Tyrrhenian samples). We were able to obtain more data on carapace length and sex at each site (N_c in Table 1) because a surplus of Nephrops were processed during the sampling activity. These data (N_c) were used to visualise variation in carapace length at each site and by sex (Fig. 1C). To statistically assess differences in size between Pomo Pits (N = 172) and other Adriatic sites (N = 135), carapace length was modelled as a function of site (a binary variable describing whether an individual was sampled in or outside of Pomo Pits) and with sex (a binary variable describing an individual as male or female).

RAD sequencing and bioinformatics

DNA extracts for each sample were sent to Floragenex (Portland, Oregon, USA) for restriction site associated DNA sequencing (RAD-seq). RAD libraries were prepared for 112 samples using the SbfI restriction enzyme and sequenced on an Illumina sequencing platform using a 2 × 100 bp approach. Raw reads were trimmed using Fastp 0.20.1 (Chen et al., 2018) and further filtered using the process_radtags program from Stacks v2.53 (Catchen et al., 2013; Rochette, Rivera-Colón & Catchen, 2019). RAD loci were built using the Stacks de novo pipeline; default parameters were used for all modules, except for -m in ustacks which was set to 3. In addition, the following parameters were added to the populations command: (i) –min-samples-overall 0.75, (ii) –min-mac 5, and (iii) –write-single-snp. Using R v4.2.0 (R Core Team, 2022), the missingno() function from poppr v2.9.4 (Kamvar, Tabima & Grünwald, 2014) was used to remove any individuals with ≥ 30% missing genotypes. Functions from dartR v2.9.7 (Gruber et al., 2018; Mijangos et al., 2022) were used to filter out loci that: (i) departed from Hardy-Weinberg Equilibrium, (ii) were in linkage disequilibrium, or (iii) were monomorphic, had a minor allele count less than five, or had all missing genotypes in a single site. Lastly, OutFLANK v0.2 (Whitlock & Lotterhos, 2015) was used to identify any outlier loci.

Genetic diversity and population structure

The gl.report.heterozygosity() function from dartR was used to calculate observed and expected heterozygosity (H_o/ H_e). The gl.report.pa() function was used to calculate the number of private alleles per sampling site.

Genetic differentiation between sampling sites was assessed by calculating pairwise values of F_st (Weir & Cockerham, 1984) using the genet.dist() function from hierfstat v0.5-11 (Goudet & Jombart, 2022). Population structure was explored using two methods: (1) a principal components analysis (PCA), and (2) a genetic clustering analysis. Prior to analysis of population structure, missing data (NAs) were imputed with the gl.impute() function from dartR using the neighbour method. The PCA was then run using the glPca() function from adegenet v2.1.3 (Jombart & Ahmed, 2011). For genetic clustering, the optimal number of genetic clusters (K) to use was determined by running the snapclust.choose.k() function and visualising the Akaike information criteria (AICc) for each K. The find.clusters() function was executed to cluster individuals into K groups based on k-means clustering. The snapclust() function was then run on the data set for the chosen K. Snapclust uses maximum-likelihood estimations to assign individual membership probabilities to each K cluster (Beugin et al., 2018). The resulting membership probabilities to each K were mean averaged per site and visualised on a projected map of the study area. Isolation-by-distance (IBD) was tested by running a Mantel test on a dissimilarity matrix of pairwise genetic (F_st) and pairwise least-cost marine geographical distances (km) using the mantel.rtest() function from ade4 v1.7-22 (Thioulouse et al., 2018).

Effective population size

The gl.LDNe() function from dartR was used to calculate effective population size (N_e), which is a wrapper around the Ne Estimator v2.1 software (Do et al., 2014). The linkage disequilibrium method was run with random mating assumed, and jackknife 95% confidence intervals were computed. Given the homogeneity of sites from the Adriatic Sea (see Results), the Adriatic sites were divided into two groups, and different minor allele frequency thresholds were tested to compare results in case of inflationary effects by rare alleles. The first group comprised all three sites from Pomo Pits (Pom1, Pom2 and Pom3), and the second group comprised of the remaining sites from the Adriatic (17I, 18II, Anc, Cgg). The rationale here was to increase the sample size in each ‘population’ and assess whether there was a difference in N_e between these two groups, since small sample sizes are likely to severely bias the estimation of N_e (England et al., 2006).

Carapace length analysis

Median carapace length in Pomo Pits was 28 mm for males (N = 95) and 24 mm for females (N = 77), while outside of Pomo Pits carapace length was 46 mm for males (N = 73) and 40 mm for females (N = 62). The model fitted to carapace length as a function of site and sex satisfied the assumptions of a linear regression (Supplementary Material S1). This model showed that, when sex is controlled for, individuals in Pomo Pits are on average 16.98 mm (41.5%) smaller than individuals outside Pomo Pits (P < 0.001). Additionally, the model showed that, when site is controlled for, males are on average 5.17 mm (12.6%) larger than females (P < 0.001).

RAD loci and genetic diversity

A total of 554 million reads were generated across 112 samples (mean = 3.78 million reads per sample). The de novo pipeline assembled 254,555 catalog RAD-tag loci, with an effective per-sample coverage mean of 4.7x (2.0x–35.7x). Filtering of individuals and loci produced a final data set of 98 individuals genotyped at 1,623 biallelic neutral SNPs; no loci were identified as outliers. At the Mediterranean sites, 0-1 private alleles were found at each site and mean observed heterozygosity was very similar across sites, ranging from 0.089−0.108. Mean unbiased expected heterozygosity values were close to observed heterozygosity, ranging from 0.094−0.113. In the single Atlantic site, the Firth of Clyde, nine private alleles were found, but heterozygosity was comparable to the Mediterranean sites.

Population genetic structure

Genetic differentiation was high between the Firth of Clyde and all Mediterranean sites, while differentiation amongst the Mediterranean sites was low or zero (Fig. 2A). The PCA revealed two distinct groups: individuals from the Firth of Clyde (Atlantic), and all other individuals from sites in the Mediterranean (Fig. 2B). This was also observed in the genetic clustering analysis, whereby K = 2 was the most likely number of ancestral populations (genetic clusters) (Fig. 2C), and this clearly showed Atlantic-Mediterranean separation into two distinct clusters (Fig. 2D). A PCA was run using only Mediterranean sites to check for any hierarchical structuring in the data (Fig. 3A). This revealed very little evidence for hierarchical structure among our sites sampled. In addition, there was little evidence for IBD across the range covered by our Mediterranean samples (Fig. 3B).

Effective population size

The N_e estimate for the Pomo Pits group was between 184 –283 depending on the minor allele frequency threshold (Table 2). For the group representing sites outside of Pomo Pits in the Adriatic, the N_e estimate was 131 –144. However, the upper confidence interval for all estimates was infinity, suggesting that, even with samples pooled to regional levels, sample sizes should be increased in future studies to ensure reliable estimates of N_e (Marandel et al., 2020).

Pomo Pits and GSA 17

The carapace length results from our study closely match the findings of Angelini et al. (2020), that is, females are on average smaller than males, and both females and males in Pomo Pits are on average much smaller than animals outside of Pomo Pits. Our genetic results indicate that Nephrops from all sites sampled in the Adriatic Sea have high gene flow between them. Likewise, using 890 samples from 27 locations genotyped at 730 SNPs, a recent Mediterranean-wide study of Nephrops did not detect any genetic differentiation within regions, including within the Adriatic (Spedicato et al., 2022: pages 493–513). This suggests that the phenotypic differences attributed to the Pomo Pits subpopulation, namely smaller mean sizes, slower growth rates and smaller mean onset at first maturity (SOM) (Angelini et al., 2020), are not explained by random genetic drift or a lack of gene flow between neighbouring Adriatic subpopulations. Moreover, Pomo Pits Nephrops have no apparent differences in genetic diversity compared to Nephrops at surrounding sites. A recent Scientific, Technical and Economic Committee for Fisheries (STECF) report showed that average Nephrops biomass from 1994-2020 is much higher in Pomo Pits (2,912) compared with Ancona (806) and GSA 18 (1,187) and outlined that Ancona and GSA 18 are at relatively lower biomass levels compared to Pomo/Jabuka Pits (Scientific Technical Economic Committee for Fisheries, 2023), both points of which accord with our genetic diversity estimations.

Our finding of regional gene flow is not unique to Nephrops across other parts of their range, nor to other closely related decapod species of the Adriatic. Using 14 neutral microsatellite markers, Pavičić et al. (2020) found the European lobster (Homarus gammarus) exhibits panmixia and comparable measures of genetic diversity in sites sampled across the Adriatic. These findings suggest that, as with Nephrops in our study, gene flow (and/or high effective population sizes) is likely sufficient to mitigate genetic drift for both species across the sites sampled. However, connectivity modelling of Nephrops has revealed the potential presence of three subpopulations in the Adriatic (Melaku Canu et al., 2021), including one along the northern Croatian coast for which no samples were available for genetic analysis in our study. This would suggest that connectivity in the region likely follows a stepping-stone model, resulting in putative genetic homogeneity across the entire Adriatic region that is caused by adequate ocean current movement, and a sufficiently long pelagic duration (1–2 months; Dickey-Collas et al., 2000) for Nephrops larval dispersal to ensure widespread admixture. Indeed, even where biophysical models strongly overestimate the realised dispersal of typical larvae, a small subset of distantly-dispersing individuals can still drive genetic homogeneity across expansive regions by inducing sufficient gene flow to offset genetic drift and effectively nullify differentiation between otherwise self-recruiting populations (Shanks, 2009; Macleod et al., 2024). Genetic connectivity at similar spatial scales has also been reported for Nephrops across Scandinavia in the north-east Atlantic, for which the authors likewise conclude that patterns of ocean currents and connectivity via larval drift are likely responsible for the observed genetic homogeneity (Westgaard, Søvik & Johansen, 2023). Although we found no genetic structuring in our sites sampled across the Adriatic, additional sampling in the eastern Adriatic will be needed to confirm panmixia of Nephrops across the entire Adriatic Sea.

In light of these findings, what could have driven the observed biological differences of Nephrops in Pomo Pits if they are not caused by genetic isolation? Other lobsters have demonstrated marked variation in SOM between discrete geographic stocks that are nevertheless highly connected genetically. Mean SOM for H. gammarus was estimated to be 18 mm CL larger for males and 31 mm CL larger for females in western Scotland than it was in eastern Scotland (Lizarraga-Cubedo et al., 2003), despite the stocks showing minimal genetic differences (Jenkins et al., 2019; Ellis et al., 2023), while the mean SOM of southern rock lobsters (Jasus edwardsii) declined markedly from 112 mm CL to just 59 mm CL along a latitudinal gradient in Tasmania (Gardner et al., 2006), despite regional genetic homogeneity (Villacorta-Rath et al., 2022). Indeed, the slower growth rates evidenced by Pomo Pits Nephrops may drive their smaller mean sizes and SOM compared to equivalent individuals in surrounding areas (Angelini et al., 2020); maturation and fecundity are often driven by age rather than size in lobsters (Gardner et al., 2006; Ellis et al., 2015), so slower growth reduces mean sizes and renders individuals smaller at a given age, reducing SOM. Environmental factors, especially temperature, are known to influence growth and SOM. Pomo Pits are considerably deeper than surrounding Adriatic areas we sampled (Table 1) and, at ∼10 °C, seawater in the Pomo depressions is typically cooler than shallower surrounding benthos (Artegiani et al., 1993). Yet, while the effect of lower temperatures on metabolism may explain the slower growth of Nephrops in Pomo Pits, this relationship is typically inversed in benthic ectotherms when it comes to SOM; for American lobsters (Homarus americanus), reduced SOM is associated with increased temperatures (Le Bris et al., 2017). Historic fishing effort may be another possible driver of lower SOM among the Nephrops of Pomo Pits. Haarr et al. (2017) modelled potential explanatory factors against spatial variation and temporal declines in the SOM of H. americanus, and, while temperature and population density were poorly correlated to SOM variations, intensity of size-based harvest selection was significantly associated, suggesting that SOM decreases were an evolutionary response to fishing pressure. Similarly, recent results from the Western Mediterranean suggest a sharp decline of Nephrops SOM which appears to be linked to increasing overexploitation of stocks (Vigo et al., 2024). However, although Nephrops stocks in the Adriatic are among the most heavily fished in the Mediterranean (Russo et al., 2019), it is not possible to assess whether spatial differences in harvest intensity may be linked to the lower SOM of Pomo Pits Nephrops; fine-scale spatial data on historic fishing effort have been lacking until recent advances in vessel tracking, and data from contemporary years suggests that other Adriatic stocks we sampled (i.e., Anc and 18II) are subject to similar levels of pressure (Russo et al., 2018). Therefore, although we cannot attribute it to any genetic divergence, any combination of (i) differences in fishing pressures over time and space, (ii) local environmental conditions, and/or (iii) intraspecific competition for space and food could have played a part in driving the exhibited variation in SOM of Pomo Pits Nephrops (Russo et al., 2018; Angelini et al., 2020; Chiarini et al., 2022; Vigo et al., 2024).

Intra-Mediterranean homogeneity

We sampled two Mediterranean sites outside of the Adriatic Sea, one in the Aegean Sea (Aeg) and one in the Tyrrhenian Sea (9I); both sites showed genetic similarity with each other and with the Adriatic sites. Our finding of minimal genetic divergence between Adriatic Nephrops and samples from adjacent sub-basins to the East and West reflects the pattern of population structure observed for the European spiny lobster (Palinurus elephas), for which samples from the Croatian Adriatic showed minimal differentiation to those from either Crete or the Balearic Sea via neutral SNPs (Ellis et al., 2023). However, H. gammarus has a similar estimated window of pelagic larval duration (PLD) to Nephrops, but does show Mediterranean sub-structuring, with Adriatic samples differentiated from those of both the Aegean Sea and the Western Mediterranean (Pavičić et al., 2020; Ellis et al., 2023). Although our depiction of homogeneity between Nephrops of the Adriatic and adjacent areas may in part reflect a paucity of samples from these sites and other regions, our analytical methods and sample sizes for Adriatic Nephrops may be sufficient to explore connectivity at this scale, since previous research has detected even subtle regional differentiation with relatively low sample sizes using 6,340 RAD-derived SNPs (e.g., Ellis et al., 2023).

Atlantic-Mediterranean divergence

At a basin-wide scale, the site sampled in the Atlantic (Firth of Clyde) was strongly differentiated from all Mediterranean sites (pairwise F_st values: 0.11 –0.13). Although only seven individuals were genotyped from this site, this pattern of Atlantic-Mediterranean divergence has also been found in previous studies of Nephrops using mitochondrial D-loop variation (Gallagher et al., 2018) and microsatellite genotypes (Gallagher et al., 2022), though homogeneity was reported using restriction fragment length polymorphisms (Stamatis et al., 2004) and allozymes (Stamatis et al., 2006), which is likely due to the lower power and resolution of these latter two molecular markers to detect genetic variation.

A pronounced phylogenetic break between Atlantic and Mediterranean populations is a characteristic common to a multitude of taxa of the north-eastern Atlantic, despite often greatly varying life history strategies (Patarnello, Volckaert & Castilho, 2007), including littoral fishes (Galarza et al., 2009), benthic echinoderms (Carreras et al., 2020) and bivalve molluscs (Wenne et al., 2022). More relevantly, other wide-ranging lobsters of the region have also demonstrated hierarchical population structure, the primary feature of which is Atlantic-Mediterranean differentiation; European lobster (H. gammarus) and spiny lobster (P. elehpas) both show strong divergence across an Atlantic-Mediterranean divide (Jenkins et al., 2019; Ellis et al., 2023). This concordant pattern amongst different lobster species using putatively neutral SNP loci suggests that, at this scale, Atlantic and Mediterranean meta-populations of all three species have been sufficiently isolated for long enough to develop differing allele frequencies and accumulate new neutral mutations through genetic drift. The drivers of this isolation are likely shaped by the topographic and bathymetric enclosure of the Mediterranean, past and present oceanographic barriers, and periodic vicariance during the Pleistocene glaciations (Patarnello, Volckaert & Castilho, 2007; Pascual et al., 2017; Jenkins et al., 2019). All of these factors serve to inhibit dispersal between the Atlantic and Mediterranean, limiting larval exchange and thus gene flow between populations in each basin (Ellis et al., 2023).