Introduction

Understanding how marine populations are connected is central to guiding conservation planning and designing effective spatial management, particularly in regions undergoing rapid environmental and coastal change such as the Mediterranean Sea (Giakoumi et al., 2013; Jahnke et al., 2017). Yet predicting gene flow in marine systems remains challenging because connectivity is shaped by the interaction between intrinsic species traits and extrinsic seascape or oceanographic processes. In species with limited adult mobility, dispersal often depends on a pelagic larval phase, whose realised contribution to gene flow is influenced by reproductive timing, pelagic larval duration, larval behaviour, self-recruitment, habitat distribution and physical transport by currents (Cowen et al., 2007; Klein et al., 2016; Ollé-Vilanova et al., 2022; Pascual et al., 2017; Sefc et al., 2020). As a result, dispersal ability does not always translate directly into genetic connectivity: some taxa with seemingly limited dispersal potential show high connectivity, whereas others with apparently greater dispersal potential exhibit pronounced population structure (Weersing & Toonen, 2009). In addition, historical demographic processes can leave persistent genomic signatures that complicate the interpretation of contemporary population structure and realised connectivity (Arranz et al., 2021; Hellberg, 2009; Selkoe et al., 2008).

The Mediterranean Sea provides an excellent setting for disentangling the processes shaping marine connectivity because contemporary oceanographic structure and historical biogeographic legacies operate together across multiple spatial and temporal scales. Present-day circulation is highly structured and seasonally dynamic, with regional gyres, jets, fronts and semi-enclosed basins that can generate dispersal barriers, asymmetric connectivity and fine-scale population structure even over short geographic distances (Cowen et al., 2007; Schroeder et al., 2012; Sefc et al., 2020). At the same time, historical events such as the Messinian salinity crisis around five million years ago and Quaternary sea-level fluctuations have shaped Mediterranean biogeographic boundaries, demographic histories and divergence times across marine taxa (Agiadi et al., 2024; Bianchi et al., 2012; Deli et al., 2019; Koblmüller et al., 2015; Kousteni et al., 2015). Thus, Mediterranean marine populations are expected to carry signatures of both contemporary dispersal processes and deeper historical events.

Comparing sympatric species offers a powerful way to test the predictability of marine connectivity. Because such species occupy the same regional seascape, many extrinsic drivers of dispersal and demographic history are shared (Dawson, 2012). Concordant patterns among species would suggest that oceanographic and biogeographic processes strongly shape realised connectivity, whereas discordant patterns would point to the importance of species-specific traits, ecological differences, or demographic stochasticity (Kato et al., 2025). Yet, although this comparative logic is central to phylogeography, it has rarely been applied at the genomic scale to closely related sympatric species. Such systems offer an additional advantage: because orthologous genomic regions can be compared within a shared evolutionary framework, they allow tests not only of whether shared seascapes or oceanography generate similar population structure, but also of whether shared environmental pressures are associated with parallel genomic landscapes of divergence, for example through selection on the same genes or pathways (Wagner et al., 2026).

The endemic Mediterranean clingfish genus Gouania (Gobiesocidae) provides an intriguing system for testing these genomic trajectories. Unlike other clingfishes, Gouania species are highly specialized for an interstitial lifestyle in intertidal gravel beaches, where adults have extremely limited mobility (Briggs, 1955; Wagner et al., 2023). Connectivity is therefore thought to depend mainly on a short pelagic larval phase (∼12 days), probably with dispersal largely constrained to nearshore environments (Macpherson & Raventos, 2006). Together with the patchy distribution of suitable gravel beach habitats, this life history should promote fine-scale population differentiation (Wagner et al., 2019, 2021).

The genus comprises five species that diversified after the Messinian salinity crisis and whose distributions broadly align with major Mediterranean ecoregions (Figure 1a) (Wagner et al., 2019, 2021). One species, G. willdenowi, is restricted to the western Mediterranean, whereas the remaining four occur in the Adriatic and eastern Mediterranean. Importantly, these include two sympatric species pairs: G. adriatica and G. pigra in the Adriatic and G. orientalis and G. hofrichteri in the eastern basin. This combination of restricted distributions, replicated sympatry, ecological similarity and low adult dispersal makes Gouania particularly suitable for asking whether shared seascapes generate predictable patterns of connectivity – and whether such predictability extends to genomic landscapes of divergence.

Figure 1

We tested this on a sympatric species pair Gouania orientalis and G. hofrichteri around Crete and the nearby island of Kythira. Crete is geographically isolated from the mainland, yet embedded in a highly structured circulation regime shaped by major jets and gyre systems and positioned at the intersection of the Ionian, Aegean and Levantine ecoregions (El-Geziry & Bryden, 2010; Spalding et al., 2007; Figure 1b, 2a). This geography creates replicated population contrasts for the sympatric species G. orientalis and G. hofrichteri around Crete and the nearby island of Kythira. In this study, we use whole-genome population data from both species to test two predictions: first, that passive larval transport should generate concordant patterns of genomic connectivity; and second, that shared environmental or ecological contrasts may produce a detectable excess of parallel genomic divergence. We then place these genomic results into a broader phylogeographic context using mitochondrial COI variation across all five Gouania species sampled throughout the Mediterranean basin.

Methods

Whole genome re-sequencing and variant calling of sympatric Gouania orientalis and G. hofrichteri around Crete

For the two eastern Mediterranean sympatric species, Gouania orientalis and G. hofrichteri, we re-sequenced whole genomes from four populations around the island of Crete and from one population of the island of Kythira, altogether leading to 46 genomes of G. orientalis and 39 genomes of G. hofrichteri (Table S1). By focusing on co-occurring species within the same oceanographic context, we aimed to disentangle shared, environmentally driven patterns of connectivity from species-specific demographic and evolutionary responses (Figure 1b).

We sequenced most individuals to a target coverage of 5X but also two individuals per population to 15X coverage (Table S1). We aligned all raw reads against the Gouania adriatica reference genome (fGouWil2.1; (Rhie et al., 2021)) using BWA v0.7.17 and extracted variants with bcftools v1.14 (Danecek et al., 2021). We filtered sites with mapping quality of zero for more than 10% of reads or root mean square mapping quality of less than 50. We further masked sites with significant (p < 0.001) strand difference and with an extremely low (<2.5 percentile) or high (>97.5 percentile) depth. We filtered heterozygous sites exhibiting significant allelic depth imbalance between reference and alternative alleles using a binomial test (PHRED score >20), affecting between 0.02% and 0.15% of heterozygous sites per individual. Further, we filtered out sites with excess heterozygosity (InbreedingCoeff < 0.2) and sites with more than 20% missing genotypes. Our final variant call sets for G. orientalis and G. hofrichteri included 5,822,506 and 3,173,859 biallelic SNPs and indels, respectively.

Lagrangian simulations around Crete

To test for correlations between oceanography and population structure for the population samples of G. orientalis and G. hofrichteri around the islands of Crete and Kythira as a null model for purely passive dispersal, we conducted Lagrangian simulations. We used the parcels package in python v3.1.1 (Delandmeter & Van Sebille, 2019; Lange & Van Sebille, 2017) based on hourly 2D current velocity data with a spatial resolution of 4.2 km from 01.11.2022 to 05.02.2025 downloaded from the Copernicus Marine Service, Mediterranean Sea Physics Analysis and Forecast (Clementi et al., 2023). Due to the coarse spatial resolution, the starting positions were not located within the current raster data. Therefore, we adapted the start position for each population for the larval dispersal modelling to the nearest raster cell centroid and used this as the new start position with longitude and latitude values (separated by comma) as follows: i) Kythira: 22.91666603088379, 36.14583206176758; ii) Petres 24.375, 35.39583206176758; iii) Mades 25.04166603088379, 35.4375; iv) Plakias 24.375, 35.14583206176758; and v) Vatos 25.58333396911621, 34.97916793823242. Throughout the whole period, we released one propagule per site each day using the fourth-order Runge-Kutta integration, following a rough estimation for the pelagic larval duration (PLD) of G. willdenowi of 13 days (Macpherson & Raventos, 2006) and recorded the location of each larva every 6 hours. Due to a lack of knowledge about the breeding season for the investigated species, we assessed seasonal drift differences between spring (March-June), summer (June-September), autumn (September-November) and winter (November-March) and we conducted a second simulation based on a PLD of 26 days.

Fine-scale genomic population structure and connectivity analyses around Crete

We then tested for gene flow and connectivity among populations calculating D-statistics in Dsuite v0.4 r42 (Malinsky et al., 2021) using the Dtrios function. For the latter, we included an individual of the other species (G. orientalis or G. hofrichteri) as an outgroup in the respective individual VCF file and adjusted the p-values for multiple testing by applying a Benjamini-Hochberg correction (Benjamini & Hochberg, 1995). In order to assess population genomic metrics, we calculated genome-wide mean estimates of nucleotide diversity (pi), Fixation index (F_ST, (Hudson et al., 1992) and mean net-between group divergence (d_A) using sliding windows of 60,000 bp with a step size of 30,000 bp. Windows containing fewer than 100 SNPs were excluded. Calculations were performed using the popgenWindows.py script from https://github.com/simonhmartin/genomics_general/tree/master and a standard calculation for d_A (following calculations and codes in (Wagner et al., 2026). We also calculated the average genome wide individual inbreeding coefficient (F_IS) using vcftools 0.1.16 (Danecek et al., 2011). To quantify the similarity of population differentiation patterns between species, we compared mean pairwise F_ST estimates for all population pairs using Spearman’s rank correlation. Correlations were calculated from the corresponding pairwise F_ST values derived independently for G. orientalis and G. hofrichteri.

For inferring phylogenetic clusters, we created a distance matrix using PLINK v1.90b6.21 (Purcell et al., 2007). Principal component analyses (PCA) were performed using PLINK2 v2.0 on genome-wide SNP datasets converted to binary PLINK format. PCA was conducted using the --pca function following filtering for minor allele frequency (MAF) thresholds of 0.01, 0.05, 0.10, 0.20 and 0.40. Eigenvectors and eigenvalues were extracted from the resulting covariance matrix and used for downstream visualization for MAF of 0.1. This led to 1,071,073 and 1,054,075 biallelic SNPs, respectively, for G. orientalis and G. hofrichteri. We investigated population genetic structure using PCA in PLINK using a MAF of 0.15 and model-based ancestry estimation with ADMIXTURE v1.3.0 (Alexander et al., 2009).

Local (parallel) genomic signals of population-specific divergence

We annotated variants from the VCF with SnpEff (Cingolani et al., 2012), using the fGouWil2.1.99 gene annotation. To improve annotation of unresolved Ensembl gene identifiers (ENSGWIG) in the fGouWil2.1.99 genome annotation, a custom annotation table was generated by combining Ensembl and NCBI gene annotations. Remaining unresolved protein-coding genes were annotated by retrieving canonical protein sequences from Ensembl and searching them against the UniProt/Swiss-Prot database using DIAMOND BLASTP (Buchfink et al., 2021). For each query, the best-scoring Swiss-Prot hit was retained and gene symbols were assigned for alignments with ≥50% query coverage and an E-value < 1 × 10⁻⁵. Existing annotations were retained and Swiss-Prot-derived annotations were only used for previously unresolved genes. This procedure provided annotations for an additional 6,554 protein-coding genes and improved the functional annotation used in downstream analyses.

First, we specifically focused on reciprocal patterns, where alleles were close to fixation in Kythira but rare or absent in mainland Crete populations, or vice versa (AF ≥ 0.95 or AF ≤ 0.05). To assess whether parallelism was also detectable at the gene level among these nearly fixed variants, we performed a permutation test in which gene sets were randomly sampled from the background set of 23,789 annotated genes while preserving the observed number of genes containing nearly fixed variants in G. orientalis and G. hofrichteri (45 and 8 genes, respectively). Significance was calculated as the proportion of permutations yielding gene overlap equal to or greater than that observed between species.

We then used two related but complementary approaches to identify genomic outlier windows. The Population Branch Statistic (PBS) (Yi et al., 2010) captures lineage-specific differentiation in the focal population while controlling for shared ancestry among populations, whereas dA identifies genomic windows with high relative divergence. To identify robust PBS outliers reflecting lineage-specific differentiation in Kythira, we calculated PBS for Kythira (PBS_Kythira) across multiple trio configurations involving Kythira and pairs of Cretan populations. For each trio-specific PBS distribution, the 95th percentile was calculated across all genomic windows. A window was classified as a PBS outlier only if its PBS value for Kythira exceeded the 95th percentile threshold in all trios, ensuring that inferred signals of branch-specific differentiation were consistent across alternative reference population combinations. Second, for each pairwise comparison between Kythira and a Cretan population, d_A outliers were defined as genomic windows falling within the 95th percentile of the d_A distribution. To account for background genetic structure within Crete, d_A windows that were also identified as outliers in any Crete-Crete comparison were excluded. Like the PBS approach, this analysis was intended to account for biases associated with population structure within the Cretan lineage, focusing on genomic regions that consistently distinguish Kythira from the broader Cretan lineage. By intersecting PBS and d_A outliers, we defined a set of high-confidence genomic windows associated with population-specific divergence in Kythira.

We then identified parallel genomic regions as those windows classified as outliers in both species under the same criteria. To assess whether overlap among candidate windows exceeded random expectations, we performed a hypergeometric test using the total number of genomic windows analysed as the background set (29,039 windows), while preserving the observed number of candidate windows in each species. Expected overlap and enrichment were calculated relative to the null expectation of random overlap. To further evaluate whether these candidate windows showed signatures consistent with selection, we calculated Tajima’s D (Tajima, 1989) within Kythira. We then tested whether candidate windows were significantly enriched for Tajima’s D values more negative than the genomic background, consistent with recent selective sweeps. To facilitate direct comparisons between populations and control for non-population specific evolutionary processes, we calculated ΔTajima’s D as the difference between Tajima’s D in Kythira and the mean Tajima’s D across Crete populations for each genomic window.

We then extracted candidate genes from these outlier windows from the Gouania adriatica annotation (fGouWil2.1.115). For each outlier window, we extracted genes overlapping a ±30 kb interval around the window midpoint on the same scaffold.

Finally, we conducted three functional enrichment analyses performed with g:Profiler (Reimand et al., 2007) using gene sets unique to (i) G. orientalis, (ii) G. hofrichteri and (iii) the overlap of windows between the two species (i.e., parallel outlier windows). This approach allowed us to investigate both independent and shared (parallel) genomic pathways underlying local adaptation. All GO enrichment analyses were performed using zebrafish GO annotations with a custom gene background (all Gouania adriatica reference genes) applying an FDR threshold of 0.05 and a Benjamini-Hochberg correction (Benjamini & Hochberg, 1995).

Mediterranean basin-wide phylogeographic reconstructions based on COI data

For assessing phylogeographic patterns across the genus, we analysed a total of 668 partial sequences of the mitochondrial cytochrome-c-oxidase subunit I gene (COI) from 23 populations, generated in previous studies (Wagner et al., 2019, 2021, 2026) and deposited on GenBank (accession numbers: MK873443–MK873539, MT299844–MT299870 and OL839338–OL944591 summarized in Table S2). We aligned and trimmed the individual COI sequences in MEGA 7.0 (Kumar et al., 2016) and applied the integrated MUSCLE (Edgar, 2004) algorithm. We visualized the geographic structure of intraspecific mitochondrial genetic diversity for each of the five Gouania species by means of a parsimony network using TCS (Clement et al., 2002) in PopART v.1.7 (Leigh & Bryant, 2015).

Demographic analysis based on whole genome data

For inferring the demographic histories of individual Gouania species, we applied the Multiple Sequentially Markovian Coalescent (MSMC2) model (Schiffels & Wang, 2020) to a deep-sequenced (∼30-fold coverage) specimen of each Gouania species (G. hofrichteri and G. orientalis from Attica, Greece; G. adriatica as well as G. pigra from Pula, Croatia and G. willdenowi from Toulon, France). The individual raw sequence files were generated in a previous study (Wagner et al., 2026) and accessed from the European Nucleotide Archive using the numbers: ERS9744577 (G. willdenowi; Toulon, France), ERS9744566 (G. pigra; Pula, Croatia), ERS9744553 (G. orientalis; Saronida, Greece) and ERS9744540 (G. hofrichteri; Chamolia, Greece). As we were lacking this data for G. adriatica, we used short-read data obtained for the creation of the reference genome (Rhie et al., 2021), by downloading a fraction of the original raw data (“fGouWil2_S1_L007_R[1&2]_001.fastq.gz”; Pula, Croatia) from GenomeArk (https://genomeark.github.io), which yielded about the same coverage as individuals from the other species. We visualized the results assuming a mutation rate of 3.5e-08 (Malinsky et al., 2018) and an estimated generation time of 2 years (own aquarium observation). For each MSMC2 run we performed 20 bootstrap replicates following the method and plotting scripts from (Lescroart et al., 2023).

Results

Lagrangian simulations predict structured passive dispersal around Crete

To generate explicit predictions for passive larval transport around the island of Crete, we performed Lagrangian passive particle drift simulations (Figure 2). Assuming a drifting duration of ∼13 days, corresponding to previous estimates of pelagic larval duration in Gouania (Macpherson & Raventos, 2006), simulations revealed pronounced spatial structure across all seasons (Figure 2b–e).

Figure 2

Kythira was largely isolated from Crete and northern and southern Cretan populations were strongly separated. Dispersal trajectories showed substantial overlap between the two northern populations, Mades and Petres, across all seasons. In contrast, connectivity between the southern populations, Vatos and Plakias, was weaker, particularly during spring (Figure 2c). Seasonal simulations further suggested limited and asymmetric exchange around Crete. During summer (Figure 2d), a small proportion of particles released from northern Crete, particularly Mades, reached the southern coastline, whereas in autumn (Figure 2e), limited northward dispersal from southern Crete, especially Vatos, was observed.

Overall, the simulations predict a spatially structured and seasonally variable dispersal regime, with strongest connectivity among northern populations, limited exchange between northern and southern Crete and weak connectivity between Crete and Kythira. The trajectories also suggest a predominantly clockwise pattern of potential passive dispersal around Crete. Using a longer drifting duration of 26 days yielded largely consistent results (Figure S1), suggesting that these predictions are not particularly sensitive to the assumed pelagic larval duration. In the following section, we evaluate this null model of dispersal driven solely by physical oceanographic processes against genome-wide population structure.

Whole-genome data reveal concordant connectivity patterns in sympatric Gouania

To test whether the passive-dispersal predictions were reflected in realised genomic connectivity, we analysed whole-genome population structure in the sympatric species Gouania orientalis and G. hofrichteri around Crete and Kythira. Across both species, we found a pronounced division between northern and southern Cretan populations, with individuals from Kythira clustering more closely with northern Crete than with southern Crete. This pattern was consistently supported by PCA, pairwise F_ST, admixture analyses and tests of excess allele sharing (Figure 3; Figures S2–S3; Table 1).

Figure 3

Table 1

Principal component analyses revealed broadly congruent population structure in both species (Figure 3a,b). In both cases, PC1 separated Kythira from the Cretan populations, whereas PC2 captured a north–south subdivision within Crete. In G. orientalis, the single individual from NNE-Mades grouped closely with NW-Petres. In G. hofrichteri, the north–south subdivision was also evident, but population clustering within Crete was more pronounced. Together, these results indicate shared geographic structure across species, characterized by strong differentiation of Kythira and a secondary north–south break within Crete.

Admixture analyses supported these patterns. In G. orientalis, Kythira was strongly differentiated from the Cretan populations (Figure S2a). In contrast, northern Cretan populations of G. hofrichteri showed admixture components also present in both Kythira and southern Crete (Figure S2b), consistent with connectivity mediated stepping-stone-wise through northern Crete. When analyses were restricted to Cretan populations only, both species showed a pronounced north–south structure (Figure S3).

Patterns of population differentiation were highly concordant between species. Pairwise F_ST estimates were strongly correlated between G. orientalis and G. hofrichteri across all population comparisons (Spearman’s ρ = 0.733, P = 0.0158), indicating broadly similar patterns of genetic structure (Figure 3c). In both species, differentiation was low among populations from the same side of Crete and generally higher between northern and southern populations. Comparisons involving Kythira showed the highest levels of differentiation. In G. hofrichteri, Kythira was less differentiated from northern than from southern Cretan populations, whereas differentiation between Kythira and Crete was more uniform in G. orientalis. Comparisons involving the NNE-Mades population of G. orientalis yielded consistently low F_ST values, but this pattern should be interpreted cautiously because only a single individual was sampled from this locality.

Tests of excess allele sharing largely mirrored the connectivity patterns predicted by passive larval drift simulations. In both species, D-statistics revealed significant excess allele sharing among northern and southern Cretan populations (Table 1), despite the inferred sister relationship between Kythira and northern Crete in the population trees. In particular, excess allele sharing between SE-Vatos and northern Cretan populations relative to SW-Plakias suggests clockwise connectivity around the eastern flank of Crete, consistent with the simulated dispersal trajectories.

Both species showed a similar relative pattern, with Kythira exhibiting significantly lower nucleotide diversity and higher inbreeding than other populations (Figure 3d; Figure S4). Individual inbreeding coefficients (F_IS) differed strongly among populations (two-way ANOVA, population effect: F_4,75 = 98.62, P = 4.37 × 10⁻²⁹), whereas the overall effect of species was not significant (F_1,75= 2.23, P = 0.140). However, a significant species × population interaction (F_4,75 = 6.72, P = 1.12 × 10⁻⁴) indicated that population-specific patterns of inbreeding differed between G. orientalis and G. hofrichteri. Genome-wide diversity metrics revealed both species-specific differences and shared population-level trends. Nucleotide diversity (π) was consistently higher in G. orientalis than in G. hofrichteri across populations (Figure S4). Despite this difference, both species exhibited remarkably similar spatial patterns of diversity. Kythira showed the lowest nucleotide diversity in both species, whereas the Cretan populations displayed comparatively higher and more homogeneous levels of diversity. The single G. orientalis individual from NE-Mades showed considerably higher nucleotide diversity than other populations, which likely contributed to the lower F_ST values observed in comparisons involving this locality.

Shared population contrasts reveal parallel and species-specific genomic divergence

The concordant population structure observed in G. orientalis and G. hofrichteri around Crete provides an opportunity to test whether shared population contrasts are associated with parallel patterns of genomic differentiation across the genome. In both species, Kythira was strongly differentiated from Crete, showed reduced nucleotide diversity and had elevated inbreeding coefficients, making it a comparable focal population across two independently evolving lineages. We therefore used Kythira-centred comparisons to ask whether divergence between Kythira and Crete involved the same genomic regions, genes, or pathways in both species, as expected if some differentiation outliers reflect local adaptation driven by shared selective pressures.

To identify genomic regions associated with Kythira-specific divergence (Figure 4; Figure S5), we quantified PBS and d_A outliers across 29,039 genomic windows. In G. orientalis, 966 windows exceeded the genome-wide 95th percentile for PBS across all three Kythira-centred trios. In G. hofrichteri, 347 windows met this criterion across six Kythira-centred trios. This difference partly reflects the exclusion of the single NE-Mades individual from the G. orientalis analysis, which reduced the number of informative Kythira-centred trios relative to G. hofrichteri. Using d_A, we detected 1,612 Kythira-specific outlier windows in G. orientalis and 2,276 in G. hofrichteri.

Figure 4

We then defined final candidate windows as those identified by both PBS and d_A. This yielded 397 candidate windows in G. orientalis (Figure 4a; Figure S5a) and 220 candidate windows in G. hofrichteri (Figure 4b; Figure S5b), whereas only 8 were shared between the two species, (i.e., parallel outlier windows). This overlap was significantly higher than expected under random overlap: 3.01 windows were expected, corresponding to a 2.66-fold enrichment (hypergeometric test, P = 0.011). These shared windows therefore represent a conservative set of parallel genomic regions associated with the Kythira–Crete contrast in both species.

At the gene level, these 397 candidate windows in G. orientalis and 220 candidate windows in G. hofrichteri corresponded to 895 and 656 unique genes, respectively (Table S3). Intersecting these sets identified 50 shared genes, which is significantly more than expected by chance (permutation test, P = 4.0 × 10⁻⁵), providing further evidence for parallel genomic divergence associated with the Kythira–Crete contrast on a gene-level.

We further examined whether parallelism was detectable at the level of nearly fixed variants. We identified 54 and 12 nearly fixed variants (SNPs and indels) associated with 52 and 14 genes in G. orientalis and G. hofrichteri, respectively (Table S3). Although no identical nearly fixed variants were shared between species, both contained Kythira-specific nearly fixed variants associated with the gene tuba (Tables S4, S5). This overlap was highly unlikely to arise by chance given the observed numbers of candidate genes (permutation test, P = 0.0334). These results therefore suggest parallelism at the gene level involving tuba, despite the absence of shared underlying variants.

We next assessed whether candidate regions showed signatures consistent with recent or localised selection by examining Tajima’s D. In contrast to G. orientalis (Figure 4c), candidate windows in G. hofrichteri were significantly shifted towards more negative Tajima’s D values relative to the genomic background (Mann–Whitney U; P = 2.91 × 10^-37), indicating a Kythera-specific excess of low frequency variants in these regions and consistent with an enrichment for population-specific selective sweeps in Kythira (Figure 4d). At the genome-wide scale, however, G. orientalis showed generally more negative Tajima’s D values across all populations, which may influence the detectability of localised selection (Figures S6, S7).

Finally, we tested whether candidate genes were enriched for specific functional categories. GO enrichment analyses of genes located in parallel windows and overlapping outlier windows revealed no significant enrichment. Likewise, species-specific candidate gene sets yielded predominantly broad and non-specific functional categories (Table S7, S8). Outlier-associated genes in both species were enriched for general molecular functions and cellular processes, including ion binding, protein binding, catalytic activity and metabolic pathways. In G. hofrichteri, enriched terms were primarily associated with phosphorylation and phosphorus-containing compound metabolism, whereas G. orientalis showed additional enrichment for protein modification, cell motility and cadherin signalling. Overall, enrichment patterns were diffuse and did not identify a single dominant biological pathway, suggesting that adaptive divergence may involve numerous genes with diverse cellular functions rather than a small number of functionally coherent pathways.

Genus-wide mtDNA variation reveals broad basin-level structure but species-specific phylogeographic patterns

Across the genus, mitochondrial COI variation revealed pronounced differences in phylogeographic structure among Gouania species (Figure 5). The western Mediterranean species, G. willdenowi and the two Adriatic species showed comparatively compact haplotype networks. G. willdenowi exhibited a shallow, star-like network dominated by a common haplotype shared among French populations, with the Messina population represented by a private but closely related haplotype (Figure 5a).

Figure 5

Similarly, G. adriatica showed a star-like network dominated by a common central haplotype (Figure 5b), whereas G. pigra showed two main haplotype groups corresponding broadly to northern and southern Adriatic populations (Figure 5c). These comparatively shallow networks are consistent with recent demographic expansions in western and Adriatic populations, possibly linked to Quaternary sea-level fluctuations and postglacial changes in shallow coastal habitat.

By contrast, the two eastern Mediterranean species, G. hofrichteri and G. orientalis, showed substantially stronger mitochondrial structure and higher haplotype diversity. However, this structure was not parallel at the level of individual populations. In G. hofrichteri, haplotypes from Crete and Kythira were largely associated with the same divergent mitochondrial cluster, separated from many mainland Greek and more northern populations (Figure 5d). In G. orientalis, mitochondrial structure was more complex and Kythira did not cluster with Crete but instead appeared closer to haplotypes from mainland/Aegean populations (Figure 5e).

Species-level demographic histories are heterogeneous across Gouania

To obtain a species-level temporal perspective on demographic patterns across Gouania, we reconstructed effective population size (Ne) trajectories with MSMC2. These analyses, based on one high-coverage individual per species, revealed substantial heterogeneity in inferred demographic histories across the genus (Figure 6).

Figure 6

Overall, MSMC2 results indicate that species-level demographic histories are only partly concordant across Gouania and cannot be readily explained by either deep phylogenetic relationships (Figure 6a), present-day sympatry, or geographic co-occurrence. Instead, demographic trajectories appear to reflect a combination of deeper historical processes and species-specific responses to environmental change, contrasting with the stronger concordance observed for contemporary genomic connectivity around Crete.

The two eastern Mediterranean species showed distinct histories (Figure 6b). Gouaniua orientalis exhibited relatively high inferred Ne through much of the late Pleistocene, followed by a decline towards more recent periods, whereas G. hofrichteri showed a lower trajectory over the same interval, with a decline between the late Pleistocene and the Last Glacial Maximum (LGM). Both species showed relatively stable Ne since the LGM. The Adriatic species, G. adriatica and G. pigra, showed declines in inferred Ne between the late Pleistocene and the LGM, followed by relative increases within the last 10 kyr, although their detailed trajectories and absolute Ne values differed markedly. The western Mediterranean species, G. willdenowi, showed a trajectory broadly similar to its sister species G. pigra, with a sharper but more variable recent increase.

Discussion

A central assumption in comparative phylogeography is that closely related and ecologically similar species exposed to similar environmental conditions should exhibit comparable patterns of connectivity and divergence (Dawson, 2012). In this context, the sympatrically occurring species of the genus Gouania provide an ideal system for testing this hypothesis across both broad (Mediterranean-wide) and fine (within-region) spatial scales. Indeed, across multiple levels of biological organization, sympatric Gouania species exhibited a remarkable degree of concordance. Contemporary genomic connectivity around Crete closely matched oceanographic predictions in both species and shared population contrasts were associated with a significant excess of parallel genomic divergence. In contrast, mitochondrial phylogeographic patterns and species-level demographic histories were considerably more heterogeneous. Together, these results suggest that shared oceanographic and biogeographic contexts can generate predictable patterns of connectivity and aspects of genomic divergence, whereas deeper historical processes remain comparatively species-specific.

Oceanography generates predictable connectivity patterns across sympatric species: insights from Crete

These patterns of predictability extend to smaller spatial scales, as demonstrated by our case study around the island of Crete. Both sympatric species, Gouania orientalis and G. hofrichteri, exhibit remarkably similar patterns of population structure that closely match the surrounding oceanographic regime. In particular, the pronounced differentiation between northern and southern Cretan populations, as well as the strong isolation of Kythira, are consistently recovered across both species and align closely with predictions from Lagrangian particle drift simulations.

This concordance between simulated passive dispersal and observed genomic connectivity indicates that oceanographic processes are likely a primary driver of gene flow in these species. Notably, this result challenges previous assumptions that Gouania, due to its short pelagic larval phase and interstitial lifestyle, might be largely restricted to nearshore dispersal (Macpherson & Raventos, 2006; Wagner et al., 2019, 2021). Instead, our results demonstrate that even limited pelagic durations (ca. 10–17 days) can be sufficient for oceanographic transport to structure population connectivity at regional scales. More generally, this is consistent with theoretical and empirical studies showing that even low levels of larval exchange can influence patterns of genetic structure in marine populations (Palumbi, 2003). At the same time, both species show elevated inbreeding coefficients in Kythira, consistent with strong local larval retention and suggesting that populations may persist as highly localized units embedded within weakly connected dispersal networks. Such nearshore retention is common in cryptobenthic fishes and is often linked to strong habitat specialization and settlement into spatially rare, patchily distributed microhabitats such as Mediterranean gravel beaches (Brandl et al., 2018; Wagner et al., 2019).

Despite these broadly parallel patterns, subtle but consistent differences between species indicate that species-specific traits still modulate realized connectivity. For example, G. hofrichteri exhibits slightly weaker population differentiation and a broader distribution range compared to G. orientalis, which may reflect differences in dispersal ability or larval ecology. Indeed, G. hofrichteri has a wide distribution range and also occurs within the Adriatic Sea, suggesting overall higher dispersal abilities (Wagner et al., 2019, 2021). This result fits observations in the sister genus Lepadogaster (Conway et al., 2019), which are also capable of long range dispersal via local oceanic features (Klein et al., 2016). However, within the same genus, differences in breeding season (and concomitant with that larval duration) between sister-species highlight the role of species-specific life-history traits in shaping realized connectivity patterns (Henriques et al., 2002; Tojeira et al., 2012; Wagner et al., 2017).

Together, our results indicate that while oceanographic structure generates highly predictable patterns of connectivity across species, the strength and fine-scale expression of these patterns remain modulated by species-specific biological traits, raising the question of whether such shared connectivity patterns also translate into parallel genomic responses.

Significant but limited signals of parallel genomic local population divergence

Similar patterns of population structure were associated with a detectable excess of parallel genomic divergence. By focusing on Kythira as a shared focal population across both species, we minimized environmental variation and compared genomic responses under broadly similar oceanographic and ecological conditions. Although most candidate regions were species-specific, overlap among both genomic windows and genes exceeded random expectations, indicating that comparisons involving the same geographic population (Kythira) can leave repeatable genomic signatures across independently evolving lineages.

Notably, one of the overlapping candidate genes was tuba, which encodes a core component of the microtubule cytoskeleton. Microtubules play fundamental roles in cell polarity, morphology and tissue organization (Larson et al., 2010), suggesting that some aspects of local divergence may involve shared cellular organizational processes across species. However, the limited number of overlapping genes and the absence of strong functional enrichment prevent more specific biological inference.

Despite this significant overlap, most candidate regions remained species-specific. Candidate regions were distributed across the genome and were predominantly located in non-coding regions, suggesting that local divergence is polygenic and likely mediated largely through regulatory variation rather than a small number of major-effect loci. This interpretation is consistent with the broad and largely unspecific functional enrichment patterns observed, as well as the limited overlap of candidate regions between species. Similar patterns have also been observed in our recent study of convergent adaptation across Gouania, where clear phenotypic convergence was accompanied by weak and diffuse genomic signals, likely reflecting the combined effects of polygenic architectures, recombination and evolutionary divergence among lineages (Wagner et al., 2026).

Differences in the extent of divergence between species further suggest an important role for demographic history. Gouania orientalis exhibited stronger population structure, higher genome-wide π, and substantially more nearly fixed variants than G. hofrichteri. Moreover, nearly fixed variants in G. orientalis were frequently clustered within genomic regions, whereas those in G. hofrichteri were fewer and more dispersed. Together with differences in genome-wide Tajima’s D, these patterns indicate contrasting demographic histories that likely influenced the accumulation of genomic differentiation, consistent with the expected effects of genetic drift on allele frequency dynamics (Kimura, 1968).

Subtle ecological differences may also contribute to the predominance of species-specific candidate regions. Although both species occur in sympatry and experience broadly similar environmental conditions, subtle differences in microhabitat use could expose populations to distinct selective pressures. Comparable small-scale habitat differentiation has been documented among Adriatic Gouania species (Wagner et al., 2023). However, the limited functional overlap between species suggests that similar ecological challenges may often be mediated through different genetic pathways.

Together, these results indicate that similar environmental and oceanographic contexts can generate repeatable genomic responses across independently evolving lineages, but that this predictability is limited. Most genomic divergence remains species-specific and appears to reflect a combination of polygenic architecture, demographic history and lineage-specific ecological responses.

Species-specific demographic histories contrast with parallel contemporary connectivity

In contrast to the strong concordance observed for contemporary connectivity and aspects of genomic divergence, basin-wide mitochondrial phylogeographic patterns and species-level demographic histories were considerably more heterogeneous. While the COI data support a broad contrast between relatively shallow western/Adriatic structure and deeper eastern Mediterranean structure, they do not show a simple one-to-one correspondence between sympatric species in the placement of individual populations. These mitochondrial patterns and the MSMC trajectories therefore provide broad historical context for our genomic analyses, but also highlight that single-locus phylogeographic structure does not necessarily mirror genome-wide connectivity at fine spatial scales.

Interestingly, the observed patterns closely align with basin-specific palaeoenvironmental histories and are most parsimoniously explained by Quaternary sea-level oscillations, with the Adriatic Sea providing a particularly clear example. During glacial lowstands (ca. 15,000 - 20,000 years ago), the sea level was around 120 m lower than today and large parts of the northern Adriatic were exposed (Lambeck et al., 2014). This drastically reduced the available shallow coastal habitats, leading to population contractions. Consistent with this, G. adriatica exhibits a pronounced star-like haplotype network indicative of recent population growth, matching the re-establishment of the modern Adriatic coastline after the last glacial maximum. In contrast, G. pigra shows a stronger north–south phylogeographic structure, suggesting either more persistent refugial structure within the basin or a more complex recolonization history. Similar postglacial (northern) expansion dynamics have been reported for other benthic taxa in the Adriatic Sea, indicating that these patterns likely reflect shared responses to basin-specific habitat dynamics rather than species-specific processes alone (Koblmüller et al., 2015; Ledoux et al., 2018; Sefc et al., 2020).

The eastern Mediterranean species further illustrate the limits of phylogeographic predictability. Both G. orientalis and G. hofrichteri exhibit elevated haplotype diversity and deeper mitochondrial structure, consistent with greater long-term habitat stability in the eastern basin. However, the placement of the populations of Crete and Kythira differed substantially between species, demonstrating that even closely related taxa occupying similar environments can retain distinct historical signatures. Beyond these basin-scale effects, intraspecific structure also aligns with major Mediterranean ecoregions (Spalding et al., 2007) and oceanographic fronts (El-Geziry & Bryden, 2010), further supporting the role of persistent seascape features in shaping connectivity. For example, the clear break between populations from Messina and the French coastline in G. willdenowi is consistent with restricted dispersal across major Mediterranean ecogeographic regions, namely the Tyrrhenian and Algero-Provencal basins (Spalding et al., 2007).

Together, the basin-wide analyses reveal a striking contrast with the contemporary genomic results around Crete. Whereas oceanographic conditions appear capable of generating repeatable patterns of connectivity and even partially parallel genomic divergence across species, mitochondrial phylogeographic structure and long-term demographic trajectories are substantially less predictable. These deeper historical signatures likely reflect species-specific responses to Pleistocene environmental change, demographic stochasticity and the persistence of historical population structure. Consequently, predictability in Gouania appears strongest for contemporary processes and progressively weaker when viewed across longer evolutionary timescales.

Conclusion

Our study demonstrates that the predictability of marine population structure varies across evolutionary timescales. At contemporary spatial scales, shared seascape features and oceanographic processes generate remarkably concordant patterns of connectivity across sympatric Gouania species. Moreover, shared population contrasts were associated with a significant excess of parallel genomic divergence, indicating that common environmental and demographic contexts can leave partially repeatable signatures across independently evolving lineages. In contrast, mitochondrial phylogeographic patterns and species-level demographic histories were considerably more heterogeneous, highlighting the increasing influence of species-specific historical processes over longer timescales.

For conservation management, these findings emphasize that efficient planning requires integrating oceanographic predictions with species-specific biological information and empirical genetic estimates of realized connectivity (Jahnke et al., 2017; Paterno et al., 2017). Our results further show that genomic data are essential for assessing the evolutionary and adaptive potential of populations at local scales, particularly in rapidly changing environments such as the Mediterranean Sea (Boulanger et al., 2022). This is especially relevant for highly specialized species such as Gouania, whose gravel beach habitats are increasingly threatened by beach nourishment and coastal development (Carević, 2020; Carević et al., 2022; Wagner et al., 2023). Accordingly, understanding patterns of connectivity and local genomic divergence will be critical for predicting demographic resilience and preserving evolutionary potential in this endemic Mediterranean genus.

Data availability

All bioinformatic analyses were implemented using workflow-based pipelines developed in Snakemake v7.25.0 (Mölder et al., 2021), and downstream analyses, visualization, and statistical analyses were conducted in Jupyter notebooks. The complete analysis workflows and scripts used in this study are publicly available at GitHub (https://github.com/maxwagn/GouaniaPhylogeographyCrete) and, for the particle-drifting simulations, at https://github.com/arianedr/Particle_Simulation/blob/main/Particle_Drifting.ipynb. Archived versions of all local analysis scripts, notebooks, intermediate files, and supplementary datasets are available through Figshare (https://doi.org/10.6084/m9.figshare.32686272). Raw whole-genome sequencing data generated for this study have been deposited in the European Nucleotide Archive (ENA) under accession ERR7917043-ERR7917126 (Table S1). Mitochondrial COI sequences analysed in this study are publicly available through GenBank under accession numbers MK873443-MK873539, MT299844-MT299870, and OL839338-OL944591 (Table S2).

Acknowledgements

The authors thank NOBIS-Austria for the financial support (NOBIS-Scholarship to M.W.). M.W. has also been supported with a DOC-Scholarship by the Austrian Academy of Science during his PhD and a YUFE4Postdoc-Scholarship, co-funded by the EU under Marie Sklodowska-Curie Grant Agreement No. 101081327. H.S. acknowledges support from the Research Fund of the University of Antwerp”.

Additional files

Supplementary tables.

Supplementary figures.

Additional information

Funding

European Union (101081327)

• Maximilian Wagner

University of Antwerp (University Research Fund)

• Hannes Svardal

NOBIS Austria (NOBIS-Scholarship)

• Maximilian Wagner

Austrian Academy of Sciences (DOC-Scholarship)

• Maximilian Wagner