Parallel connectivity and genomic divergence but heterogeneous demographic history in sympatric Mediterranean clingfishes

Reviewer #1 (Public review):

Summary:

This article taps into the very interesting phylogeographic situation of two sympatric species of clingfishes sharing the same distribution and environment around Crete and the island of Cythera. Basically, it shows that the population structures of these fishes are influenced in a parallel manner by seascape, low dispersal, and potentially drift or selection despite a different phylogeographic history. This parallelism is looked after in a detailed manner at the genome level. I am very enthusiastic about this extremely well-constructed and very cleverly designed study of this natural "common garden" evolutionary duplicate situation, something sufficiently rare to be underlined.

The authors have produced complete genomes for the two species and their five population samples, from which they are able to conduct up-to-date data analyses. The text is clearly written, without too much jargon, and the options chosen for the various analyses and bioinformatic pipelines are sufficiently detailed so it is rather easy to follow what they've precisely done, something which is alas not that frequent in comparable studies.

Strengths:

The structural part of their study is really very convincing, with results showing that despite very small differences, the population structures of the two species conform quite well with what could be inferred from larval dispersal modeled according to passive particle drift. The parallelism is striking, despite minute differences, and despite quite different population sizes for the two species.

Weaknesses:

After that, the authors tried to identify a set of outlier loci whose distribution doesn't conform to the main population differentiation. This search for outliers is made according to classical methods based on Fst or its derivates, like the program PBS which compares the Fst values in trios on sliding windows along the genome. The study is well conducted, namely taking into account false positives and false discovery rates in a conservative manner.

Assuming that some environmental variables differ between the five locations where they have samples, the authors hypothesized from this that similar environmental pressures give similar patterns, potentially affecting pro parte similar places in the genome of the two species. However, there is a blind spot in their analysis inasmuch as it seems that recombination and linkage are not taken into account. It is well known that recombination rates are very variable along the chromosomes, going from recombination hot spots to stretches of very low levels of crossovers. It is well known as well that a variety of phenomena like background, purifying, and sweep selection are quite sensitive to the recombination rates, this having important bearings on local variation of a series of variables like Fst, D, Pi... . Hence, their conclusions about a direct action of environmental pressures may largely be challenged, especially when it is to look for functionality of tightly linked genes, and should be taken as the last hypothesis to be retained when the others can be ruled out.

There are probably sufficient levels of conservation and synteny in teleost fish and a sufficient number of species where recombination maps exist so that the authors can reconstruct a map for their species, at least partially, permitting to use methods explicitly taking recombination rate variation along the genome into consideration (like for instance DILS: Demographic inferences with linked selection 2021 Mol Ecol Res) to challenge their adaptationist conclusions, all the more given the fact that the question of the eventual nature of differential environmental pressures cannot be addressed with extant data.

(3) Moreover, there is now a considerable amount of literature which deals with what is coined "islands of differentiation" or "islands of speciation" or "barriers loci". These genomic segments coincide very often with low recombination regions, and some of them are quite often shared in multiple pairs of closely related species. It seems that their experimental set up (several closely related species) is ideal to easily derive the landscape of genomic architecture of divergence in the genus, this permitting to see where the conserved islands of differentiation stand, how much they are conserved or not in the different species, and how this matches or not with their PBS peaks, and by the way, allowing a more direct comparison with similar landscapes published in other species. They have everything at hand, and this will be a very valuable addition to this article that could hence become a more fascinating paper.

Reviewer #2 (Public review):

Summary:

The authors analyzed genetic population structure of two clingfish species in the eastern Mediterranean Sea. They used a genome-wide DNA sequence dataset representing several populations of both species and found general patterns of isolation that agree with likely patterns of dispersal and described parallel signatures suggesting genomic adaptation. The results support that similar species follow the same evolutionary trajectories when they evolve in the same ecological and geographical context and when meta -populations are subject to the same constraints with respect to dispersal.

Strengths:

The findings as such include a population genetic analysis of taxa for which data are lacking. Conclusions on population subdivision as well as genomic divergence suggestive of genomic adaptation are well supported. The analysis of the genetic data is according to established standards and useful.

Weaknesses:

The idea to apply drift simulations to study dispersal in pelagic fish larvae makes sense; I do however, wonder to what extent the neutral drift scenario applies to the species under study, as an earlier paper by the authors ( J Fish Biol 2020 Nov 3; 98(1): 64-88) emphasizes a strong near-coastal retention of larvae. Accordingly, the drift simulations are in agreement with the observed population structure, but I note that a simpler isolation-by-distance model and dispersal primarily along coastlines would equally agree with the data. A demonstration that drift simulations contribute to a refined model of population subdivision would require a denser sampling of populations and measures of drift between these. Ideally, this should include a demonstration that areas where drift is reduced coincide with genetic discontinuity in the absence of deeper areas of the sea. This is not to say that the approach as such is not interesting, but the conclusions towards this goal are not very well supported by data.

One aspect in the study should be clarified: between-species gene flow is common between closely related species of fish. If this occurs, this would dramatically affect the interpretation of parallelism and convergent evolution at the genomic level. Hybridization should be ruled out or discussed. The two focal taxa are not sister taxa according to a phylogenetic tree (Figure 6) based on only a few genetic loci. Even if a more complete dataset is missing for most relevant taxa to recalculate this phylogeny, I would like to see an assessment of the degree of separation between the two taxa at a genome-wide level. The authors mention that the genetic variants that are subject to parallel evolution are not identical. This is promising and should be highlighted more in case hybridization is likely.

Reviewer #3 (Public review):

Summary:

This manuscript asks whether biogeographic patterns are similarly predictable across temporal scales in two sympatric clingfishes (Gouania). By integrating oceanographic dispersal simulations, population genomics, mtDNA, and demographic inference, the authors find broadly concordant contemporary connectivity between species but less consistent deeper phylogeographic and demographic histories. The study addresses an interesting question, although several aspects of the genomic analyses and their interpretation need further attention.

Strengths:

The manuscript has a strong comparative design and an interesting biological system. Comparing sympatric species provides a useful framework for asking how consistently shared environmental processes translate into similar evolutionary patterns. The integration of oceanographic, population-genomic, phylogeographic, and demographic evidence across temporal scales is a particular strength. The most convincing result is that contemporary connectivity appears more predictable across species than deeper phylogeographic and demographic history.

Weaknesses:

The main new dataset is whole-genome resequencing of 85 individuals (46 G. orientalis, 39 G. hofrichteri) from five localities on Crete and Kythira. These reads were mapped to the chromosome-level Gouania willdenowi assembly (fGouWil2; publicly labeled as G. willdenowi but referred to throughout as the G. adriatica reference) produced by the Vertebrate Genomes Project (Rhie et al. 2021). No new reference genome was generated here, which limits the genomic contribution. Ideally, each focal species would have its own reference.

Mapping both species to a third congener also creates an important reference-bias issue. Figure 6a places G. orientalis closer to the reference lineage than G. hofrichteri, so the more distant species is expected to map less well and yield fewer callable sites. That is the direction observed: G. hofrichteri yields nearly half as many called variants (3.17 M vs 5.82 M). This could reflect real diversity, differential mapping, or both. Per-species mapping rates, missingness, and callable-site counts should therefore be reported, and outlier comparisons repeated using sites callable in both species.

Coverage is also low and uneven: 69 individuals were targeted at 5X and 16 at 15X. This is probably adequate for broad population structure, but more problematic for heterozygosity- and frequency-based statistics, as hard calls at ~5X can bias heterozygosity, nucleotide diversity, FIS, and nearly fixed variants. The single 15X G. orientalis from NNE-Mades, which has unusually high diversity and affects FST, illustrates the issue. Key results should therefore be checked using genotype likelihoods, common-depth downsampling, or the higher-coverage subset. Table S1 should also report realized depth and missingness and clarify the discrepancy with the Methods, which state that two individuals per population were sequenced at 15X. Agreement among PCA, admixture, and FST is not independent reassurance because all use the same hard-called genotype matrix (see Lou et al. 2021; Helmkampf et al. 2025).

Much of the remaining dataset comes from previous studies or public resources: the COI sequences, four of five MSMC2 genomes, the VGP genome, and Copernicus currents. Thus, the main new contribution is the Crete-Kythira population-genomic dataset. The strongest result is that contemporary connectivity is broadly concordant between species, whereas deeper phylogeographic and demographic histories are less predictable. The geographic mismatch among datasets, however, limits direct comparisons across scales.

The least secure result is parallel genomic divergence, which is emphasized in the title. The null model for parallel divergence is problematic. The enrichment test treats 29,039 windows as independent even though 60-kb windows overlap by 30 kb. Adjacent windows are therefore not independent, and both species also share the same reference architecture. The significance of 8 shared windows among 397 and 220 outliers may consequently be inflated. This should be tested with a null model that preserves genomic autocorrelation, such as non-overlapping blocks or permutations of larger genomic segments. Otherwise, the parallelism claim should be softened, particularly in the title.

The gene-level analysis raises the same concern and also contains several inconsistencies. The background is 23,789 genes in the Methods, but 21,274 in Table S3; the Methods specify {plus minus}30 kb around window midpoints whereas Table S6 appears to use 120 kb; and the nearly fixed-variant analysis uses 45 and 8 genes in the Methods versus 52 and 14 in the Results/Table S3. These numbers need to be reconciled. I would also avoid placing much weight on a single shared tuba gene without an appropriate null model.

The MSMC2 analysis uses one diploid genome per species but is interpreted as species-level demographic history. Because these genomes come from different localities and the manuscript itself shows substantial within-species structure, species and locality are confounded. These results should be framed as exploratory histories of the sampled genomes/populations unless additional individuals are analyzed.

The link between oceanography and genomic connectivity is also mainly visual. Because this is a central aim, it should be tested quantitatively, for example by comparing a Lagrangian connectivity matrix with genetic distance while accounting for geographic distance. Similarly, the correlation between pairwise FST values across species should use a matrix-based permutation test rather than a standard Spearman p-value based on 10 non-independent comparisons. The Methods explicitly state that correlations between oceanography and population structure are tested, but no such analysis is reported.

A few additional points need attention. The NNE-Mades G. orientalis is excluded from some analyses but retained in others despite strongly affecting results. The particle-drift model assumes passive 2D advection and omits larval behavior and vertical movement, which should be stated. Mitochondrial-nuclear discordance could also reflect stochasticity, introgression, or sex-biased dispersal. Finally, several inconsistencies need correction, including duplicate Figure S6/missing S7, four versus five populations in the Figure 2 caption, the ENA accession range, and several typos.