Highly restricted dispersal in habitat-forming seaweed may impede natural recovery of disturbed populations

Cystoseira sensu lato (Class Phaeophyceae, Order Fucales, Family Sargassaceae) forests play a central role in marine Mediterranean ecosystems. Over the last decades, Cystoseira s.l. suffered from a severe loss as a result of multiple anthropogenic stressors. In particular, Gongolaria barbata has faced multiple human-induced threats, and, despite its ecological importance in structuring rocky communities and hosting a large number of species, the natural recovery of G. barbata depleted populations is uncertain. Here, we used nine microsatellite loci specifically developed for G. barbata to assess the genetic diversity of this species and its genetic connectivity among fifteen sites located in the Ionian, the Adriatic and the Black Seas. In line with strong and significant heterozygosity deficiencies across loci, likely explained by Wahlund effect, high genetic structure was observed among the three seas (ENA corrected FST = 0.355, IC = [0.283, 0.440]), with an estimated dispersal distance per generation smaller than 600 m, both in the Adriatic and Black Sea. This strong genetic structure likely results from restricted gene flow driven by geographic distances and limited dispersal abilities, along with genetic drift within isolated populations. The presence of genetically disconnected populations at small spatial scales (< 10 km) has important implications for the identification of relevant conservation and management measures for G. barbata: each population should be considered as separated evolutionary units with dedicated conservation efforts.


Results
A total of 503 individuals from 15 sampled populations ( Fig. 1; Table 1) was amplified using the set of twelve microsatellite markers newly developed; 383 individuals successfully amplified at six markers at least, with up to two alleles per marker. Out of these 383 genotyped individuals, eight pairs of individuals (2.1%) were presumed to be clones. These were observed across several distinct sites (two pairs in site 5, in site 9 and in site 12; one pair in site 7 and in site 8). In addition, these pairs of identical individuals were labeled with sequential-or very close-numbers. As we can not rule out that the same individual was sampled twice, we decided that one unique genotype per clone was retained for further analysis. This resulted in a final dataset composed of 375 unique genotypes using twelve microsatellite markers.
Quality control of the newly developed markers. The quality and variability of the microsatellite markers were tested for population genetic studies of G. barbata based on a sampling comprising 15 sites in three Seas, the Ionian, the Adriatic and the Black Seas. All loci were found to be polymorphic (5-39 alleles, Table 2), with H s estimates varying from 0.242 (Cysbar2) to 0.736 (Cysbar5). F IS values ranged from − 0.710 (Cysbar6) to 0.501 (Cysbar3), with significant deviations from Hardy-Weinberg Equilibrium observed for seven loci (Cysbar3-Cysbar6, Cysbar9-Cysbar11). Missing data were present in all subsamples and in ten out of the 12 microsatellite markers (Supporting Information SI1), likely attributed to the presence of null alleles.
Of the twelve microsatellite markers tested, MicroChecker suggested the presence of null alleles in one to five subsamples (out of 15) for seven loci (Cysbar1 to Cysbar5, Cysbar10 and 11). However, using unilateral exact binomial tests, heterozygosity deficiencies were likely explained by the presence of null alleles in only Cys-bar3 (p-value = 0.133). MicroChecker did not detect any stuttering or short allele dominance (SAD), however, the regression approaches detected significant SAD in Cysbar5 (ρ FIS = − 0.530, p-value = 0.0002; ρ FST = − 0.591, p-value = 0.00003). Stuttering was not detected following De Meeûs et al. 67 approach. The analysis of 375 G. barbata individuals genotyped with twelve microsatellite markers showed significant linkage disequilibrium between two loci, cysbar6 and cysbar12, after Benjamini and Yekutieli adjustment. We removed three problematic microsatellite markers, Cysbar3 and Cysbar5 for which technical issues were observed, and Cysbar6, in linkage disequilibrium with Cysbar12, from the final dataset.
Once removing these three markers (N = 375 individuals, nine microsatellite markers, 6.7% missing data, SI1), we analyzed the stability of F IS and F ST from one locus to the other following the determination key proposed by De Meeûs 68 . A strong and significant heterozygote deficiency was observed with f = 0.066 (p-value = 10 -3 ). The F IS variation across loci was moderate (Fig. 2), but a high F ST variation was observed among loci, the Pearson correlation between F IS and F ST was not significant (r(F IS ,F ST) = 0.257 p-value = 0.397) and the standard error of F IS (StrdErrFIS = 0.022) was two times lower than for F ST (StrdErrFST = 0.046). A strong and significant correlation  Table 1. The map was created using the R package mapdata 105-107 . Strong genetic structure at large spatial scale. High (Fig. 3a), in line with the strong and significant heterozygote deficiencies observed, partly explained by the presence of null allele, but also likely explained by a Wahlund effect, which justified a clustering approach. This result is congruent with the Principal Component Analysis along the first axis (PCA; SI2) and the genetic distances ( Table 3). Note that different values of K were explored (from two to 18) and K = 2 provided the best meaningful result, not only using the Evanno et al. 69 method (ΔK K=2 = 4558.9, ΔK K>2 < 1.5), but also by crosschecking these results with PCA and F ST estimates. Further increases in K were provided in SI2.
These two lineages, i.e. Mediterranean samples and Black Sea samples, required being analyzed separately. . The F IS variation across loci was higher than the variation observed among all subsamples (Fig. 2), the Pearson correlation between F IS and F ST was not significant (r(F IS ,F ST) = − 0.133 p-value = 0.733), and the standard error of F IS (StrdErrFIS = 0.066) was 1.5 higher than for F ST (StrdErrFST = 0.043). No correlation was observed between the number of missing data and F IS values (ρ = − 0.017, p-value = 0.517). Further, within the Ionian and Adriatic Seas, a strong structuration was revealed by the clustering analyses. Again, different values of K were explored (from two to ten, results provided in SI3) and K = 3 (Fig. 3b) provided the best meaningful results (ΔK K=3 = 1551.4, ΔK K=4 = 112.6, ΔK K=2, 5-10 < 1, but also by crosschecking these results with PCA and F ST estimates). Geographic groups comprised (i) sites 1-2, (ii) sites 3-4 and (iii) sites 5-6. A similar genetic structure was detected by PCAs, provided in SI3. Table 2. Characterization of twelve microsatellite markers developed in Gongolaria barbata (15 sites sampled, N = 375 unique genotypes). N a : Number of alleles, H S : Nei's genetic diversity within subsample; F IS : fixation index and its p-value from U test for heterozygote deficiency or excess (**p-value < 0.01, ***p-value < 0.001, non-significant for all the other loci). The three loci bolded were removed from the final dataset.  13   www.nature.com/scientificreports/ showed mixed ancestry to three genetic backgrounds, suggesting some local gene flow. In addition to mixed ancestries observed in sites 9-12, some individuals were assigned to one of the two genetic backgrounds with an ancestry rate higher than 0.8. Further increases in K did not result in new meaningful clusters (ΔK K < 11, except ΔK K=2 = 101.3, results provided in SI4). Clustering from PCA (provided in SI4) appeared similar to clustering revealed by Structure.
Within sea estimates of dispersal. Significant correlations were observed between genetic differentiation and oceanographic distances among sites, either using the shortest path, i.e. the shortest distance between two sites, or the longest path, i.e. following the coastlines, in the Black Sea, but only using the shortest path in the Adriatic Sea (Fig. 4). In the Adriatic, we observed a slope b = 0.121, which, once translated, suggested a relatively small neighborhood size N b = 8 individuals and an immigration from neighbors within the same sea of N e m = 1.32 individual per generation and subpopulation. With an averaged N e = 78.55, and based on an estimated number of populations between 10 and 100 in the Adriatic Sea, the estimates of effective population density ranged from 0.006 to 0.06 individuals/km 2 , i.e. very sparse populations, and the dispersal (δ) was estimated between 186 and 587 m, i.e. very restricted dispersal. In the Black Sea, based on the longest path that better explained the regressions, the slope b = 0.0238 was translated into a neighborhood size N b = 42 individuals and an immigration from neighbors within the same sea of N e m = 6.69 individuals per generation and subpopulation. With an averaged N e = 50.76, and based on an estimated number of populations between 100 and 1000 in the Black Sea, the estimates of effective population density ranged between 0.014 and 0.140 individuals/km 2 , i.e. very sparse populations, though higher than in the Adriatic, and the dispersal (δ) was estimated between 129 and 408 m, i.e. very restricted dispersal as estimated in the Adriatic.

Discussion
In the present study, we investigated the genetic diversity and structure of fifteen remnant populations of G. barbata sampled in three seas, the Ionian, the Adriatic and the Black Sea, aiming to evaluate whether the natural recovery of disturbed populations may occur without the help of restoration actions. The use of microsatellites for conservation genetic purposes has undoubtedly be a major advance, providing several advantages over allozymes as nuclear markers 65,70 . But a variety of potential pitfalls need to be carefully considered when using microsatellites, especially regarding null alleles 67,68,71,72 . The determination key 68 applied on the microsatellite markers specifically developed for G. barbata allowed us interpreting data with heterozygote deficits and discriminating technical to demographic causes. Facing G. barbata population decline, low level of intra-population genetic variation would be expected due to recent demographic bottlenecks in wild populations. Interestingly, the polymorphism observed in G. barbata is on the same range than that observed not only in other Cystoseira species 54,57,61 , but also, in other Fucales [73][74][75] . Across the fifteen G. barbata sampled sites, we found significant heterozygote deficiency, which is common in Cystoseira 54,56,58 . However, in most previous published studies, markers with potential technical issues were not removed nor discussed, excluding the possibility of direct comparisons among species regarding the possible causes explaining these deviations.
Here, the technical issues of the microsatellite markers were ruled out and allele frequencies corrected to consider the presence of null alleles, leaving the origin of the significant deviations from HW expectations to be    55 . In the present study, however, departures from HWE in G. barbata cannot be explained solely by the presence of inbred individuals: inbreeding was detected in a single site located in the Adriatic Sea. In addition, we did not detect selfing from the analyses of the multilocus heterozygosity distributions. In Cystoseira s.l., fertilization is oogamous, with thousands of biflagellate male gametes meeting female gametes in the vicinity of hermaphroditic conceptacles grouped within terminal receptacles. This reproduction mode strongly suggests that selfing may occur in these species if gametes are self-compatible. But, we could not find solid evidence of self-compatibility in the literature, or direct citation of an experimental study specifically testing for it 55,56,80 , a common assumption that should therefore deserve further investigation. Given this, mating among relatives would be more likely, without solid evidence of self-compatibility in Cystoseira s.l., facilitated by the very low dispersal of the fertilized eggs 29 , leading to clustering of related individuals 55 . This said, deviations from HW expectations may also be attributed to a Wahlund effect: given the very limited dispersal abilities of the species 29,30 , mixing of populations that differ in allele frequencies would result in significant heterozygote deficiencies in the pooled population, a scenario that is very likely, especially if more related individuals settle close to each other so that populations are made of a sum of different families. One of the most expected results was certainly the strong genetic structure observed across the three seas, with the convergence of all methods towards nearly no gene flow between the Black Sea and the Ionian and Adriatic Seas. Strongly differentiated populations have been observed not only in other Cystoseira s.l. species overall their distribution ranges 54 , but also in other canopy-forming seaweeds [81][82][83] . The genetic distinction of the Black Sea populations is very likely explained by a lack of exchange of migrants through the Bosphorus strait, consistent with previous finding based on mtDNA markers 62 . A very similar structuration pattern was found over the same sampled sites in Tritia neritea, a gastropod characterized by a direct cycle without a pelagic larval stage 84 . The restricted exchange of water masses and the difference in water density salinity and temperature between the water masses of the Mediterranean and the Black Seas across the Bosphorus strait has already been reported as a main driver of genetic structure in various species [84][85][86][87][88][89] . Besides the environmental gradients, transition zones are often zones of secondary contact of two lineages, i.e. a zone of contact between lineages previously isolated into glacial refugia 90 . Current data make it difficult to establish the historical demographic contexts of the observed genetic divergence, the identification of any barriers to gene flow that promote/maintain genetic isolation, as well as the exact spatial segregation of the two main genetic clusters in G. barbata. Given the observed differentiation between the Mediterranean and Black Seas populations, both using nuclear and mitochondrial markers 62 , and the inter-basins context, this pattern is likely resulting from divergence without gene flow. The extent to which www.nature.com/scientificreports/ these two genetic backgrounds are reproductively isolated deserves further attention which could be addressed, for instance, with a study covering the entire range distribution of the species, and especially covering the contact zones to identify barrier to gene flow 91 . The presence of reproductively isolated lineages would prevent the success of genetic rescue if the populations crossed were found too divergent to reproduce. Genome-wide data will likely help to delineate the presence of reproductively isolated lineages and test the hypothesis of one or several evolutionary species 92 . Similar and high genetic diversity was nevertheless observed in these two genetic backgrounds, supporting long-term regional persistence of large populations rather than recent expansions and isolation by distance, which are often associated with loss of diversity in fucoids 93 .
In addition to the strong differentiation between the Mediterranean and Black Seas observed in G. barbata, genetic analyses revealed high genetic structure at fine spatial scale, indicative of limited population connectivity, even between sampled sites separated by few kms. Such levels of genetic differentiation are in the range of that observed among populations of Cystoseira s.l. studied over the similar spatial scales (i.e. within seas) using microsatellite markers 11,55,56,79 , suggesting similar dispersal abilities of G. barbata to other Cystoseira s.l. The estimated dispersal distance per generation in G. barbata was extremely small, likely ranging from 129 to 587 m across the Adriatic and the Black Sea. Although these estimated dispersal distances should be carefully interpreted 94,95 , these estimations were nonetheless comparable to distances estimated from spatial autocorrelation analyses in E. amentacea 55 , and in agreement with estimated distances from field observations 29 . These estimated dispersal distance were likely smaller than the distance between two consecutive suitable habitat patches (shallow sheltered coasts and lagoons), enabling the migration to reverse the negative effect of genetic drift. Our results suggest that, so far, G. barbata isolated populations did not seem to suffer from inbreeding depression, and displayed intrapopulation genetic variation similar to other Fucales; genetic and demographic monitoring should be pursued to address the most appropriate conservation issues, e.g. assisted gene flow, minimizing non-genetic threats, or in case of small inbred populations, genetic rescue with highly divergent populations 50 .
Pairwise genetic differentiation among sampled sites were significantly correlated with geographic distances in the Black Sea and Adriatic. Such isolation-by-distance (IBD) is common in macroalgae, whatever their algal division, habitat and life cycle 96 , and usually observed in Fucales [97][98][99] . However, in E. amentacea, genetic structure was not significantly related to geographic distances; this relationship was improved when dispersal probability across available rocky habitat was estimated by simulating stepping-stone directional transport mediated by ocean currents 55,58 . Considering the very low dispersal life history of E. amentacea, dispersal in this species would be facilitated by possible transport in floating rafts. Indeed, rafting is a relatively frequent and documented mechanism of marine dispersal allowing fragments of thalli to be transported entangled in floating rafts of other algae 100 . In G. barbata, relative long dispersal could also occur through drifting of fertile individuals or ramification, bearing fertile receptacles, thanks to their aerocysts, helping in preserving inter-population genetic variation. Greater connectivity was, for instance, observed in C. compressa compared to E. selaginoides and E. amentacea, which was ascribed to the presence of aerocysts, providing a buoyancy property that E. selaginoides and E. amantacea do not possess 54 . Dispersal as fertilized zygotes has also been observed in Cystoseira s.l. under strong hydrodynamic events 58,100 but would not allow such long-distance dispersal. Such properties would enhance connectivity among suitable habitat patches in G. barbata, especially within the Black Sea as compared to the Adriatic Sea. Indeed, while differences in the range of geographic distances among sampled populations may partially explain this result, this finding could also reflect: (1) the presence of more free-living individuals within the Black Sea as compared to the Adriatic, (2) the higher homogeneity of water masses in the Black Sea; the oceanographic circulation is composed of a counter-rotating (cyclonic) Rim Current surrounded by many predominantly anticyclonic coastal eddies 101 while the three cyclonic gyres (the north, the central and the southern Adriatic sub-gyres) and the physical and biogeochemical features operating in the Adriatic Sea creates natural barriers for dispersion 102 and (3) larger effective population size in the Black Sea, in line with higher abundances in the Black Sea than in the Mediterranean Sea, these three hypothesis being non-exclusive. Nevertheless, despite potential occasional long-distance gene flow through rafting or drifting, strong genetic differentiation among populations together with homozygotes deficiencies suggest very restricted connectivity among natural populations of G. barbata, a result coherent among Cystoseira s.l. studies [54][55][56]58 . G. barbata populations are likely connected via neighboring subpopulations in the form of stepping stones dispersal, as supported by significant IBD. Considering the limited dispersal estimates found here, our results call into question the natural recovery of disturbed populations if suitable patches of habitat are more than 1 km far apart.
The restricted connectivity observed among G. barbata subpopulations suggests that each population represents a different evolutionary significant unit, giving the highest conservation priority to this threatened species 103 . Conservation should be designed to preserve the distinctive genetic diversity here observed. However, conservation efforts must be complemented by a reflection considering all the management units and their exchanges through gene flow, e.g. genetic rescue by assisted gene flow, to counteract genetic diversity erosion within each evolutionary significant unit. Other active strategies should be considered to restore the connectivity between these currently isolated populations, fostering the natural recovery of degraded populations (e.g. limiting the regression causes) or actively restoring lost populations 15,17,19,20 . Regarding the very restricted dispersal abilities, more genetic data on Cystoseria based on a sampling not only covering the large but also the very fine-scale distribution ranges and across time are needed to conserve the iconic Mediterranean and Black Seas marine forests and upscale reasoned ecological restoration programs.

Methods
Sampling collection and DNA extraction. Cystoseira sensu lato corresponds to three monophyletic clades recognized as three different genera 62,104 . The two clades outside the Cystoseira sensu stricto clade were recently better redefined, so that the formerly-known Cystoseira barbata was first renamed Treptancha barbata 104
Microsatellite development. Genomic DNA of 15 individuals from different locations in the Mediterranean and Black Seas was sent to the Ecogenics GmBH company (http:// www. ecoge nics. ch), commissioned to identify and develop a set of twelve Single Sequence Repeats (SSRs), i.e. microsatellite markers. Briefly, size selected fragments from genomic DNA were enriched for SSR content by using magnetic streptavidin beads and biotin-labeled CT and GT repeat oligonucleotides. The SSR-enriched library was analyzed on a Roche 454 platform using the GS FLX Titanium reagents. A total of 24,770 reads were sequenced with an average length of 474 base pairs. Of these, 3008 reads (12.2%) contained a microsatellite insert with a tetra-or a trinucleotide of at least six repeat units, or a dinucleotide of at least ten repeat units. Suitable primer design was possible in 1997 reads. Twelve markers (0.6%) showed a correct amplification and variability by polymerase chain reactions (PCR).
Following PCR optimization and polymorphism analysis of each of these twelve primer pairs, a multiplex procedure was performed using Multiplex Manager 108 and four combinations of multiplex were tested. Based on banding patterns checked on agarose and polyacrylamide sequencing gels (2% and 8%, respectively), two multiplex reactions of six microsatellites with forward primers labeled with ABI fluorescent dyes ( Table 2)  Testing the dataset quality. Before further analysis, we checked for clonal propagation across all samples by estimating the probability of identity (PI, i.e. the probability that two independently sampled individuals within a mating population share an identical multilocus genotype by chance) using GenAlEx v6.502 109 . With a PI < 1.8 10 -5 per population, any identical genotypes found in the dataset were presumed to be clones. These clones could either belong to a unique genet or be two representatives of the same clone that have spread. To discriminate between these alternative hypotheses, we carefully checked the sampling site(s) of the two individuals presumed as clones and their label numbering, used as a proxy of the distance between these individuals.
For each pair of loci, linkage disequilibrium (LD) was tested using the G-based test implemented in Fstat 2.9.4 110 , which was proven to be the most powerful approach to combine tests over all subsamples 111 . The 66 tests conducted (12loci × (12-1)/2) were not independent and p-values were corrected using the Benjamini and Yekuteli procedure 112 to obtain FDR (false discovery rate) corrected individual p-values.
MicroChecker software 113 was used to detect the presence of null alleles and eventual scoring errors (due to stuttering and/or large allele dropout/short allele dominance). We reported the null allele frequencies estimated using the Bookfield's second method provided by MicroChecker, and computed (1) null alleles frequencies expected under Hardy-Weinberg Equilibrium, and (2) the number of null alleles expected based on these frequencies and the number of individuals genotyped. We then tested the adjustment between observed and expected numbers of missing data 67,71 with a unilateral exact binomial (H1: there is less missing data observed than expected). Additionally, we tested for genotyping errors due to the occurrence of short allele dominance (SAD) using unilateral (ρ < 0) Spearman's rank correlation tests between allele size and F IS. These tests were also performed with F IT , aiming to increase the power of statistical tests 71 . Genotyping errors due to the occurrence of stuttering were tested following De Meeûs et al. 67 .
Analyses of genetic diversity. Allele frequencies, the average number of alleles (N all ), allelic richness (A r ), gene diversity (H S ) and the fixation index (F IS ) were estimated for each study site using Fstat 2.9.4 110 . Briefly, A r is the expected number of alleles in a random subsample of size g drawn from the population 114,115 . This rarefaction method allows evaluation of the expected number of different alleles among equal sized samples drawn from several different populations. Genetic diversities between two groups were compared using bilateral tests with 10,000 randomizations implemented in Fstat 2.9.4 110 . We tested for heterozygote deficiency across loci and subsamples using the U test 112 implemented in Genepop 4.7.5 116 .
To test for possible inbreeding occurring within each sample, individual Multi-Locus Heterozygosity (MLH) was calculated for each individual. In a population with variance in inbreeding, inbred individuals are less heterozygous (i.e. lower MLH) than non-inbred ones. Inbreeding variance generates Identity Disequilibria (ID)-i.e. correlations in homozygosity across loci, a measure of departure from random associations between loci 117 . IDs were measured for each sampling site by calculating the g 2 parameter and its standard deviation using the RMES software 118 . The g 2 parameter measures the excess of double heterozygotes at two loci relative to the expectation under a random association, standardized by average heterozygosity, providing a measure of genetic association www.nature.com/scientificreports/ and inbreeding variance in the population, free of technical biases. To test for the significance of g 2 , random re-assortments of single-locus heterozygosities among individuals were tested using 1000 iterations. Selfing rate within each subsample was also estimated using RMES 118 . To test whether the selfing rate significantly differed from 0, both the unconstrained and constrained likelihood options were run. We then computed Δdev = 2 * (ln l(unconstrained) − lnl(constrained)) and compared it to a chi-square with one degree of freedom 118 .
Analyses of genetic structure. First, we estimated population genetic structure through the approach developed in De Meeûs 68 based on the Wright's fixation indices 119 , F IS , which measures inbreeding of individuals relative to inbreeding of subsamples, F ST , which measures inbreeding of subsamples relative to total inbreeding, and F IT , which measures inbreeding of individuals relative to total inbreeding. In other words, F IS measures the deviations form Hardy-Weinberg equilibrium (i.e. local panmixia), F ST measures genetic differentiation between subsamples, and F IT reflects the combination of both. Briefly, this approach is based on the examination of the F-statistics stability across loci and subsamples to interpret data with heterozygosity deficiencies. Heterozygosity deficiencies can arise from null alleles, stuttering, short allele dominance, allele dropouts (i.e. technical problems), Wahlund effect (the juxtaposition of several groups with different allele frequencies), and/or selfing and sib-mating. By examining the (1) missing data in the dataset (across loci and subsamples), (2) 120 , i.e. f, θ and F for F IS, F ST and F IT estimators, respectively. We computed the 95% confidence intervals of F-statistics using the StrdErrFIS and StrdErrFST computed by jackknife over subsamples, for each locus, and we used the 95% confidence intervals computed by bootstrap over loci. Correlations between F IS and F ST were measured using Pearson's correlation coefficients in R 107 . All these analyses (parameters estimates, jackknife, bootstrap) were performed using Fstat 2.9.4 110 and R 107 .
We further tested a correlation between F IS and the number of missing data 71 : the significance of correlations was tested with a unilateral (ρ > 0) Spearman's rank correlation test in R 107 .
In addition to the aforementioned approach, genetic structure among sampling sites was depicted using more classical tools, aiming to look whether geographic differentiation exists. We used three different approaches: (1) F ST estimates (2) Principal Component Analyses (PCA) computed on the matrix of genotypes, and (3) an individual-based Bayesian clustering method. Comparing results from analyses using different statistical approaches allows us to make solid assumptions about our data: it allowed us to cross check the outputs of a model-based approach with strong priors and hypotheses (HW equilibrium, no linkage between markers) as implemented in Structure 2.3.4 121,122 , with methods based either on distance with few (nearly no) assumptions on data (i.e. PCA) or on distance matrices among sets of individuals that fall within predefined population samples. Estimates of F ST were calculated using FreeNA with the ENA correction 72 . To test for genetic isolation among subsamples, we used the G-based test over all loci. The PCAs were carried out using the R package adegenet 1.4-2 123 . The individual-based Bayesian clustering analysis was performed with the software Structure 2.3.4 121,122 . For each value of K (ranging from 1 to 18), 30 replicate chains of 200,000 Markov Chain Monte Carlo iterations were run after discarding 25,000 burn-in iterations. An admixture model with correlated allele frequencies was applied with a priori information on sample origin. For each value of K, 30 output files were produced by Structure 2.3.4 with slightly different ancestry proportions; using Clumpp 124 , we determined individual ancestry proportions (q-values; i.e. assignment probabilities) that best matched across all replicate runs for each K value. We visualized individuals' assignment in the R software 107 .
We estimated and tested hierarchical F-statistics using the HierFstat package 125 with three levels, the Mediterranean vs. the Black Seas subsamples (level 1), the Ionian, Adriatic and Black Seas (level 2), and all subsamples (level 3).
Estimates of larval dispersal. Finally, we tested for the presence of genetic Isolation-By-Distance (IBD) by investigating the relationship between using Cavalli-Sforza and Edward's chord distance 126 (Table 3) and the natural logarithm of oceanographic distances (Ln(Dgeo) 127 . Pairwise oceanographic distance between sampled sites was estimated within each sea as the shortest oceanographic distance according to the 'Vincenty (ellipsoid)' method using the geosphere package 128 . However, as the shortest path does not reflect the sea hydrodynamics 129,130 , pairwise oceanographic distances between sampled sites were also estimated following the coastlines, calculated by summing all the straight lines connecting each cape on a map at 1:500,000 on QGis 3.18. For the Black Sea, which is circular, having the choice between following the movement against the clock or clockwise to calculate the coastline distance between two populations, only the shortest distance was kept.
Regressions between genetic and geographic distances as described above were investigated within the Adriatic on the one hand and within the Black Sea (excluding subsample 11) on the other hand, using Fstat 110 . With significant IBDs, the IBD slope allows estimating the neighborhood N b = 1/b and the number of immigrants from neighbors N e m = 1/(2πb) = 4πD e σ 2 , with D e the effective population density and σ 2 the average of squared axial distances between adults and their parents 131,132 . To estimate D e and σ 2 , we first estimated the effective population size N e within each Sea (Adriatic and Black Seas): we measured the non-random association of alleles at different loci within F1 population 133 using a bias correction 134,135 implemented in NeEstimator v2 software 136 . The P crit was adjusted following recommendations 135 and random mating was applied. Estimating N e allows computing effective population densities D e , as the averaged N e multiplied by the number of subsamples within a sea divided by the surface of the sea. The number of populations in the Adriatic Sea was estimated as follow: (1) Gongolaria barbata stands for 10% of the Cystoseira s.l. observed in the Mediterranean Sea 32 , (2) in the Adriatic, www.nature.com/scientificreports/ ca. 100 Cystoseira s.l. populations were recorded, (3) we estimated that 10% would be Gongolaria barbata, so that the number of populations in the Adriatic would be close to 10. Despite this estimation, we thought very low, we decided to test different number of populations, ranging from 10 to 100. In the Black Sea, the number of populations was estimated to be 100 63 , but we also tested different number of populations, ranging from 10 to 100. Once D e computed, a rough proxy of dispersal (δ), in order of magnitude, can be computed as δ ≈ 2 (1/ (4πDe × b)) −1/2 (i.e. N e m = 4πDeσ2 < = > δ ≈ 2 (1/(4πDe × b)) -1/2126 ).