A total of 51 mtIGS haplotypes (GenBank: JN084346-96) were identified in the 352 individuals of F. ceranoides belonging to the 22 “core” populations. The median-joining network revealed three mtIGS lineages displaying nearly disjunct geographic distributions (Table 1 and Figure 1). The Bayesian phylogenetic analyses failed to recover the temporal sequence of lineage splitting, but their monophyly was supported by high branch posterior probabilities (Figure 2a). Each phylogroup was defined by one interior and widespread haplotype. Phylogroup A, composed by haplotypes A1 and 30 related ones, was present from VIG to CAM (Western sector, W), and further south in VIA. Phylogroup B, composed by B1 and 10 related haplotypes, was distributed from ANL to CED (North-Western sector, NW), although a few B haplotypes were also detected in ORT and BAR. Finally, phylogroup C, composed by C1 and 8 related haplotypes, was exclusively found from ORT eastwards to NAV (Northern sector, N), and further east in the populations of VIL, SAN and BAY. Several peripheral populations were geographically closer to populations across the phylogeographic breaks than they were to their nearest population within the same sector. For instance, ANL (NW sector) is geographically closer to CAM (W sector; ~34 km) than to RCO (NW sector, ~70 km), and the distance between CED (NW sector) and ORT (N sector; ~36,5 km) is smaller than between CED and FER (~58 km; both NW sector).

Mean allelic richness (A), Nei’s gene diversity (H _E), observed heterozygosity (H _o) and multi-locus inbreeding coefficient (F _IS,) were estimated from for the microsatellite data-set; Haplotypic richness (N _hap), haplotypic diversity (H _hap) and nucleotide diversity (π _hap) are based on the mtIGS data-set. The mtIGS lineages and haplotypes are listed for each population (coded as in Figure 1). Absolute frequencies of haplotypes are in parenthesis (if N > 1). *GenBank accessions of the private haplotypes of the populations of VIA, VIL, SAN and BAY 15 .

Globally, only the three dominant haplotypes (A1, B1 and C1) plus two derived ones (A2, B2) were shared among at least two populations. The remaining 46 haplotypes were population-specific and among these, 12 represented non-singleton variants. Many W and NW populations harboured private haplotypes in relatively high frequencies. This pattern was apparent even in populations located inside the same drainage systems, such as UMI and ULL (Ria de Arousa, W sector), or RCO, BET, ARE and FER (Artabro Gulf, NW sector). H _hap was high in the W (H _hap = 0.717) and NW (H _hap = 0.671) sectors due to the presence of most of the local haplotype radiations, but π _hap was considerably higher in the former. Contrastingly, H _hap was very low in the N sector (H _hap = 0.166), but further east the populations of VIL and SAN possessed private, C1-derived haplotypes in relatively high frequencies.

The results of the AMOVAs showed that the sectors considered accounted for about 83% of the molecular variance of NW Iberian F. ceranoides (Table 2). Within the W, NW and N sectors, 32%, 63% and 2% of the molecular variance was accounted for by the molecular differences among respective populations. The mismatch distributions did not reject the spatial expansion of phylogroups A (P = 0.957) and C (P = 0.750), but failed to support the expansion of phylogroup B (P = 0.009; Figure 3).

P values are based on 1000 permutations.

*P > 0.05.

The seven microsatellite loci revealed a total of 76 alleles in the 352 “core” individuals genotyped (6–23 per locus), although 42% of these had global frequencies below 0.01. H _E was rather variable among populations, ranging from 0.144 (ULL; W sector) to 0.563 (BAR, N Sector) (Table 1). Approximately half of the populations exhibited significant heterozygote deficiencies. Among sectors, H _E and H _hap were not correlated (Additional file 1: Figure S1a and S1b). The W sector had the lowest (yet most variable) H _E (0.457) and the highest H _hap (0.717), whereas the N sector had the highest H _E (0.553) despite very low H _hap (0.166). The NW sector showed intermediate levels of diversity for both markers. FST ranged from 0.021 (VIG vs. UMI) to 0.685 (PON vs. CED), whereas D _est ranged from <0.001 (VIG vs. PON) to 0.879 (PON vs. ANL) (Additional file 2: Table S1). Within sectors, pairwise differentiation of populations was of the same order of magnitude, but more variable within the W sector (Additional file 1: Figure S1c).

The microsatellite genotypic clusters recovered with the FCA showed a remarkable correspondence with the mtIGS phylogroups (Figure 2b). The most obvious exceptions were the populations of LIR and CAM, both belonging to the W sector, whose genotypes appeared intermediate between W and NW and W and N sectors, respectively. Excluding these admixed populations, the W populations formed the most differentiated cluster among the three, as in the phylogenetic tree (Figure 2a). The STRUCTURE analyses showed a similar picture (Figure 4). Based on the ΔK ad-hoc criterion 55 the highest hierarchical level of genetic sub-division of F. ceranoides occurred between the W sector and the NW and N sectors (K = 2; Additional file 3: Figure S2). Further subdivision of genotypes into W, NW and N sectors (K = 3) represented a weaker, but nevertheless significant, level of population subdivision. Again, the individuals of LIR and CAM showed variable degree of admixture between the W and the NW and N sectors, respectively. Within these 3 major groups some sub-structuring was also evident (Figure 4). Iberian F. ceranoides could be subdivided into a maximum of 12 (stable) genetic clusters (K = 12), corresponding to smaller, less resolved geographic regions. Some admixture (or mixed-ancestry) between neighbouring clusters was pervasive, but also apparent in a few well separated population pairs (e.g. LIR & RCO, ANL & CED).

A significant IBD pattern was detected in the whole NW Iberian region (P = 0.001), as well as in the W (P = 0.009) and N (P = 0.009) sectors (Figure 5). In the W, however, the relationship was lost when the admixed populations of LIR and CAM were removed from the analysis (P = 0.135; data not shown).

The historical isolation, divergence and enduring integrity of three distinct and disjunct genetic pools at such a narrow spatial scale are as remarkable as they are puzzling. The congruence between the mtDNA matrilines and nuclear background shows that genetic differentiation in NW Iberian F. ceranoides is genome-wide, which excludes stochastic or selective sweeps as the drivers for the drastic mtDNA shifts. Also, the depth of genomic differentiation between phylogroups is high and characteristic of vicariant sub-populations that have long been diverging independently through accumulation of de novo mutations, drift and lineage sorting. Currently there are no large, estuarine-free regions along this shoreline and therefore distance per se cannot be invoked as a major contemporary factor accounting for the much deeper genomic differentiation among than within NW Iberian genetic sectors. The mean distance separating any two neighbouring populations (~32 km) and neighbouring populations across phylogeographic breaks (<40 km) are of the same order of magnitude, and several peripheral populations are geographically closer to populations across the phylogeographic breaks than they are to their nearest population within the same sector.

The parapatric divergence of Iberian F. ceranoides seems rather unlikely. The divergence between interior mtIGS haplotypes (A1, B1 and C1; 4–5 mutations, K2P = 0.010–0.014) is about one fifth that between F. ceranoides and F. vesiculosus (21–23 mutations, K2P = 0.043–0.046), which are estimated to have diverged between 0.73 and 3.77 million years ago ( 56 , estimate based on a calibrated phylogeny of 13 nuclear genes). Assuming a constant molecular clock, the divergence of phylogroups of F. ceranoides could date back to 165–867 ky BP. This indirect estimate is almost certainly inaccurate, but implies that phylogroups started to diverge well before the last glacial maximum. The high levels of genetic endemism and diversity of the Iberian region agree with the expectations for long-term persistence in glacial refugia. However, the global effects of past climate and sea level changes also affected Iberian refugial areas 57 58 59 . The periodic transgressions and regressions associated with the expansion/melting of land-based ice-sheets caused continuous geographical rearrangement of near-shore habitats 60 61 , estuaries included. The geographic locations of NW Iberian estuaries in the past were different from today, as were to some extent the climatologic, oceanographic and hydrologic regimes 58 59 . Even discounting significant changes in the density and location of NW Iberian estuaries throughout past millennia, it appears rather improbable that any oceanographic feature acting as a powerful demographic filter could have remained relatively static in approximately its current (and very narrow) positions during such a long and dynamic period.

The contemporary genetic sectors in Iberian F. ceranoides are more likely to result from the regional reassembly of vicariant phylogroups into their current distributions following a period of independent, mostly allopatric divergence in past contracted areas of occurrence. This scenario implies the historical fragmentation and divergence of Iberian F. ceranoides in separate refugia (refugia within refugia), the subsequent expansion of these vicariant phylogroups along contiguous shorelines, and very limited gene-flow across meanwhile established secondary contact zones. The mismatch analyses and the distribution of the mtIGS variation are compatible with this scenario of independent range expansions. Within each sector, the only widespread haplotypes are the interior (presumably ancestral) haplotypes, whereas derived haplotypes (presumably younger) are typically restricted to single populations. In the terrestrial realm, similar expansion/contraction cycles have been invoked to explain the mostly parapatric distribution of lineages and sister-species in a range of Iberian taxa currently displaying relatively continuous distributions (reviewed in 62 63 64 ).

The regular climatic changes and the dynamic shoreline/drainage geography across the Pleistocene glacial/interglacial cycles 57 58 59 have probably played an important role in the regional range dynamics of Iberian F. ceranoides, but it is impossible with the present data to establish the specific drivers and its spatio-temporal contexts with detail. Our data nevertheless permit us to speculate on the temporal sequence of colonization of NW Iberia if haplotypic diversity within each sector is assumed to represent a good proxy for the time since colonization. The arrival of phylogroup C to Northern Galicia is presumed to post-date the establishment of phylogroup A in Western Iberia. Indeed, most populations of phylogroup A harbour private haplotypes in relatively high frequencies and some even local haplotype radiations. Remarkably, the number of haplotypes found in these 8 estuaries sampled along a coastline sector as small as 150 km far exceeds the number found in central and northern Europe (N_pop = 12; 15 ). The current distribution of this phylogroup probably represents a stable interglacial rear-edge that may have experienced southward expansions during colder periods such as the last glaciation. In contrast, the entire N sector is dominated by the ancestral haplotype C1. Derived C1 haplotypes were found further East, in the populations of VIL, SAN and BAY, and may indicate that phylogroup C arrived more recently to NW Iberia from an eastern Cantabrian refugium. Finally, the 6 populations forming relict phylogroup B are probably close to their refugial distribution centred on the Artabro Gulf. In the absence of complementary demographic information, however, the chronology of colonization (leading to secondary contacts) remains highly speculative.

In this contraction/expansion scenario, the enduring integrity of the fine-scale phylogeographic structure within F. ceranoides can only be explained by very limited gene-flow across phylogeographic breaks. Migration may be particularly depressed there due to the presence of contemporary oceanographic barriers to dispersal, or be as low as elsewhere and simply reflect the inherent low vagility of the species.

The most distinctive feature of NW Iberia coastline is its “rias”, drowned river valleys formed during the marine transgression that followed the last glaciation. These rias are generally divided, based on their orientation, size and main geomorphological elements, into “Lower Rias” (between VIG and XAL), “Middle Rias” (between LIR and CED), and “Higher Rias” (eastwards of ORT), with the transitions at Cape Fisterra and Cape Ortegal, respectively. The distribution of mtDNA lineages of F. ceranoides matches these subdivisions remarkably well, although establishing a link between geographic and genetic subdivisions remains difficult, as their circumscriptions are not based on distinctive climatic, hydrological or oceanographical features that could be relevant in terms of species ecology or dispersal.

The movement of buoyant, surface-drifting seaweed rafts is constrained by near-shore circulation patterns, wind, coastline morphology, tidal currents and river plumes, although the prevailing shelf/slope circulation patterns also play a role in offshore transport. It is highly challenging at the study scale to track the movements and fate of reproductive drifters leaving/arriving estuaries and therefore to directly estimate migration rates between populations within and across sectors. However, several lines of evidence suggest the absence of any specific seascape feature generating persistent (year-round) physical discontinuities matching the location of the observed genetic breaks. In NW Iberia, circulation patterns are complex and seasonally variable 65 66 . During the Autumn–Winter downwelling season, SW winds prevail and a poleward current flows over the Western slope, with inter-annual variability in intensity and penetration into the Cantabrian Sea. During the Spring–Summer upwelling season, prevailing winds shift to become predominantly NE/N oriented, and an east/southward current develops over the shelf. These characteristic patterns are intermittently dominated by short-scale meteorological events that regionally intensify or reverse circulation during short periods in each season 67 . The physical continuity of this coastline is well illustrated by the fate of the oil spilled by the Prestige tanker 250 km west off Cape Fisterra. The leaked (buoyant) fuel reached Cantabria (830 km from the sinking point) in just 17 days and spread along the Spanish shoreline from Vigo to the Basque Country 68 , i.e., throughout and beyond the region studied here.

The available genetic evidence seems to confirm this. In the presence of barriers, genetic discontinuities would be expected to be concordant in location and eventually in depth, between co-occurring species sharing similar dispersal characteristics. However, most population genetic studies of shallow coastal biota focus on organisms with planktonic dispersive stages, and/or have very poor sampling resolution in the studied region. Demes from both the Western and Northern coasts of Galicia have been analysed in mussels 69 , stalked barnacles 70 , spider and swimming crabs 71 72 , flatfish 73 , octopus 74 and direct-developing cephalopods 75 . Their common characteristic is the absence of genetic structure, indicating widespread connectivity by marine currents over NW Iberia. The single exception is the low dispersal, ovoviviparous snail Littorina saxatilis, which showed mild subdivision north and south of Cape Fisterra 76 . The authors propose that this break coincides with a putative ecological barrier but do not exclude the alternative hypothesis that it represents a secondary contact zone where allopatrically diverged sub-populations are being homogenized very slowly. The swift spread of the invasive seaweed Sargassum muticum throughout Galicia, Cantabria and Portugal, shortly after its first detection in Asturias (1980) and in Galicia (1986; 77 ) also fails to support oceanographic barriers to drifting seaweed dispersal in this area.

The absence of prominent seascape or ecological barriers matching the phylogeographic breaks implies two things. First, that their positions may be rather contingent and simply reflect the idiosyncratic sequence of (re)colonization of NW Iberian estuaries by the three phylogroups during past range expansions; and second, that the unusually sharp genetic discontinuities at secondary contact zones are maintained despite normal migration rates.

A pattern of approximate stepping stone expansions originating in different refugia could promote the formation of genetic sectors even in the absence of dispersal barriers. This possibly reflects a relatively short viability of reproductive structures after frond dislodgement, or a density-dependent effect. F. ceranoides is dioecious and therefore effective estuarine colonization requires at least one male and one female fertile frond to be in close contact after dispersal (while synchronously releasing gametes) to produce in situ the foundational zygotes that will eventually initiate a new population. Entangled mats of drifting viable male and female F. ceranoides are more likely to form near established populations that may export significant amounts of freshly dislodged drifters. Successful colonization across intermediate and larger distances surely occurs, as exemplified by the haplotype sharing of ANL and CED, or the colonization of Norway (across the North Sea) and Iceland, but it is probably much less frequent. Anyway, the process will self-reinforce: if populations in the interior of a sector go extinct, favoured recolonization from nearby sources will preserve the pre-existing genetic pattern.

Rare effective inter-estuarine dispersal, while allowing the colonization of unoccupied estuaries, implies that gene-flow has little effect in counteracting differentiation between fully established populations, slowing or preventing any progress towards migration-drift equilibrium. The remarkable genetic homogeneity of northern Europe, which was colonized post-glacially, clearly demonstrates the lack of gene-flow from the interior of the species range, where F. ceranoides exhibits considerably more diversity 15 52 . This study goes further in demonstrating that the effects of gene-flow are remarkably insignificant even at the shortest possible scale—between consecutive estuaries. Indeed, most populations in the W and NW sectors harbour private haplotypes in relatively high frequencies, including those located inside the same drainage systems. If migration is typically so low that even consecutive populations within sectors are genetically independent, it will also be unable to readily homogenize pre-existing differentiation remaining from the colonization process.

The primacy of colonization history over ongoing gene-flow reflects a poorly connected metapopulation system regulated by dispersal processes that are only effective within very restricted spatial and temporal windows. The apparent paradox of extremely limited population connectivity (here reported within and between sectors) despite an evident colonization potential (at least in the long-term) suggests extreme density-barrier effects. While the initial colonizers may reproduce relatively free of competition and contribute disproportionately to the genetic make-up of establishing populations, once populations become fully established the large disparity in the number of residents (descending from the colonists) and subsequent immigrants act as a demographic buffer against changes in allele frequencies 46 51 .

Organisms such as F. ceranoides that inhabit patchy habitats and possess the capacity for rapid population growth and habitat saturation (compared to immigration rates) may be particularly prone to density-barrier effects. Fucus spp. are fecund, fertilization success is typically very high and recruitment in the vicinity of parental plants can be very efficient 78 79 . Monopolization of local space by marine Fucus species can be fast compared to their spread along unoccupied discontinuous shores 80 . Compared to its marine congeners, F. ceranoides may occupy vacant estuaries even faster and to a larger extent. The enclosed and sheltered nature of its habitat should improve the number and success of spawning events 81 . Furthermore, F. ceranoides is a structural species that frequently forms monospecific (and often compact) belts within its particular tidal/salinity range. Contrary to the open shore where saturated communities compete for space, its estuarine habitat is free of similar competitors, potentially contributing to increased growth rates, densities and space monopolization. The sheltered and densely occupied habitat of F. ceranoides together with the short range of gamete dispersal of Fucus spp. can also explain the observed heterozygote deficiencies as reflecting local inbreeding. Other possible causes for such deviations are unlikely, as small sample sizes appear unimportant when comparing with previously analysed larger samples 52 , no recurrent amplification failure (null alleles) was observed in any locus or population, and no evidence supports the co-existence of distinct but sympatric sub-populations within estuarine patches (Wahlund effects).

The limited lineage admixture reported here indeed suggests that F. ceranoides has relatively short temporal windows of opportunity between estuarine colonization and saturation during which rare immigration can potentially result in detectable gene-flow. Incipient reproductive isolation (pre- or post-zygotic) can add to density effects and further depress gene flow between divergent phylogroups (e.g. 41 43 ). The admixed nuclear background of LIR and CAM show that these divergent lineages of F. ceranoides can interbreed and should not be regarded as incipient species, although hybridization is common between Fucus species 82 83 . Importantly, the shallow population structure within sectors shows that gene-flow is similarly reduced between more closely related populations. Unrecognized biophysical, ecological or reproductive barriers remain valid (and mutually non-exclusive) alternatives to pure demographic effects, but are probably not the most important factor underlying the apparent lack of gene-flow across sectors that maintains the parapatric structure of Iberian F. ceranoides.

In revealing that populations of Fucus ceranoides are to a large extent bounded by extinctions and (re)colonisations, our results challenge earlier assertions that marine species are mainly unified by gene-flow. In this seaweed, connectivity estimates based on allele frequency divergence (e.g. FST and related measures) are inflated (within sectors) or depressed (across secondary contact zones) to a great extent by historical colonization processes. Despite its evident parapatric structure, an IBD correlation was unexpectedly recovered in NW Iberian F. ceranoides, showing that IBD patterns can arise where a causal relationship between distance and gene-flow is missing. Such spurious correlations have been noted in other highly structured species 84 85 , and confirm that gene-flow may not be the decisive factor underlying many significant associations between geographic and genetic distances.

The “core” populations of F. ceranoides used in this study were collected in the estuaries of all major rivers between Vigo (VIG, SW Galicia) and Navia (NAV, W Asturias), in NW Iberia (N = 22; Table 1; Figure 1). This corresponded approximately to an array of discrete but neighbouring populations with an average proximity of about 33 (±17) km. Four additional Iberian populations that are not contiguous to this “core” population set were included in some analyses; VIA (northern Portugal), VIL (eastern Asturias), SAN (Cantabria) and BAY (southern France). These and also 3 “core” populations—NOI, RCO and POR—were previously analysed by Neiva et al. 15 52 . All collection sites typically contained monospecific belts of F. ceranoides attached to hard substrata and were exposed to steep salinity fluctuations throughout the tidal cycle. At each site, 5–10 cm tips of apical vegetative tissue was excised from 16 individuals sampled along a 100–200 m linear transect or random walk; tissue samples were individually stored dehydrated in silica-gel crystals until DNA extraction. To keep sample sizes constant, a random subsample of 16 individuals was used from the previously analysed populations.

Genomic DNA was extracted from approximately 10 mg dried tissue using the Nucleospin® Multi-96 plant kit (Macherey-Nagel Duren, Germany), according to the manufacturer’s protocol. Individuals were sequenced for the mitochondrial 23 S/trnK intergenic spacer (mtIGS, 15 ), and genotyped for 7 microsatellite loci developed for congeners 87 88 89 that had shown polymorphism in Iberian F. ceranoides 52 . Primer sequences and amplification details were the same as in Neiva et al. 15 52 . Amplified fragments were run in an ABI PRISM 3130xl automated capillary sequencer (Applied Biosystems, CCMAR Portugal). MtDNA sequences were aligned, proofread and edited in GENEIOUS 3.8 90 . Microsatellite alleles were manually scored in STRAND 91 using the 350 ROX™ size standard (Applied Biosystems).

The geographic distribution of the mtDNA variation was mapped and the genealogic relationships of haplotypes were inferred using the median-joining algorithm implemented in Network 4.5 92 . A phylogenetic tree for the mtDNA sequences was reconstructed with MrBayes 93 using the best-fit model of nucleotide substitution and using Fucus vesiculosus as the outgroup (GenBank no GQ385125). Among the 88 models evaluated in jModeltest 94 95 , the HKY + G model was selected based on the Akaike information ranking. Two parallel Metropolis-coupled Markov chain Monte Carlo searches, each with four chains, were run for 2 × 10⁶ generations, sampling every 100 generations. The number of substitution rates (Nst = 2) and among-site rate variation (Rates = Gamma) were set according to the substitution model selected, leaving the remaining options as default. 10⁵ generations (1000 trees) were discarded as burn-in, and the remaining 38000 used to produce 50% majority-rule consensus trees and to calculate branch posterior probabilities.

Nucleotide (π _hap) and haplotypic (H _hap) diversity within populations and inferred mtIGS phylogroups (see RESULTS) were calculated with DNASP 5.10 96 . Summary statistics of the microsatellite genetic diversity, including microsatellite allele frequencies, mean allelic richness (A), Nei’s gene diversity (H _E), observed heterozygosity (H _O) and inbreeding coefficients (F _IS), were calculated with GENETIX 4.05 97 . The partitioning of genetic variation between and among the mtIGS sectors was examined with molecular analyses of variance (AMOVA) in ARLEQUIN 3.1 98 . The significance (P > 0.05) of the fixation indices was calculated after 1000 permutations of individuals within sectors. For each phylogroup, the occurrence of recent spatial expansions (assuming constant deme size) was tested with ARLEQUIN 3.1 98 , fitting the implemented model to the observed mismatch distribution. Significance was assessed with 1000 permutations.

The microsatellite population structure was assessed with both individual (genotype based) and population (allele-frequency based) approaches. First, the degree of congruence between the mtIGS structure/phylogeny and its nuclear background was visually inspected with the factorial correspondence analysis (FCA) implemented in GENETIX 4.05 97 . Population genetic structure was further examined with a Bayesian, model-based genetic admixture analysis implemented in STRUCTURE 2.3 99 100 . Individuals were combined into one dataset for analysis, without any a priori population assignments and admixture was allowed. Each number of assumed populations (K, set sequentially from 1 to 14) was ran five times using a burn-in of 200000 iterations and a run-length of 1000000 iterations, which was determined to be sufficient to have consistent results. The “true” number of K was inferred both from the posterior probability of the data, hereafter referred to as L(K), and following the ΔK choice criterion of Evanno et al. 55 , better suited to detect heterogeneous patterns of dispersal or co-ancestry.

Pairwise F _ST (θ; 101 ) was estimated with GENETIX 4.05 97 and pairwise D (D _est; 102 ) was estimated with SMOGD 1.25 103 . Isolation by distance (IBD) was evaluated for full and sub data-sets using reduced major axis regressions of pairwise estimates of population’s genetic differentiation against minimum marine distance, as measured in Google Earth 5.1. The statistical significance of the genetic and geographic associations was assessed with Mantel tests (1000 randomizations, P < 0.05) in IBDWS 104 .