Populations of F. distichus were sampled on both sides of the Pacific and the Arctic/Atlantic, including the understudied California rear-edge, the trans-Pacific Aleutian Arc, and the Greenland coast. At each site, tissue samples (16 < n < 24) were, whenever possible, collected every other meter along linear transects. Multiple collections were obtained at specific sites (Maine, California, Alaska) where locally distinct eco-morphotypes and/or discontinuous tidal heights were apparent. Genomic DNA was extracted from silica-dried tissue using the Nucleospin 96 Plant II kit (Macherey-Nagel, Germany). Twenty microsatellite loci previously developed for other Fucus species (Engel et al., 2003; Wallace et al., 2004; Perrin et al., 2007; Coyer et al., 2009) were amplified in a geographically diverse panel of populations to assess cross-amplification success and levels of polymorphism (data not shown). The seven most polymorphic microsatellites were selected and used to generate multi-locus genotypes for all individuals of F. distichus (see Supplementary Table S1 for details). The 23S/trnK mitochondrial intergenic spacer (mtIGS, a marker with good phylogeographic resolution in Fucus; Neiva et al., 2010, 2014; Coyer et al., 2011a; Laughinghouse et al., 2015), was screened in the same populations. Amplified fragments were run in an ABI PRISM 3130xl automated capillary sequencer (Applied Biosystems) at CCMAR, Portugal. Microsatellite alleles were manually scored in STRAND (Veterinary Genetics Laboratory, University of California, Davis; http://www.vgl.ucdavis.edu/STRand) using the 500 LIZ size standard (Applied Biosystems). Sequences were aligned and edited with GENEIOUS 4.8 (Biomatters; http://www.geneious.com).

For microsatellite data, summary statistics of genetic diversity (allele frequency, mean allelic richness [A], unbiased gene diversity [H _E], observed heterozygosity [H _O] and inbreeding coefficients [F _IS]) were calculated with GENETIX 4.05 (Laboratoire Génome, Populations, Interactions, Université de Montpellier II; http://kimura.univ-montp2.fr/genetix) and the R package standArich v1.0 (available at http://alberto-lab.blogspot.pt/p/code.html). Pairwise population differentiation (Jost’s D, Weir and Cockerham’s Θ) was calculated with the R package DiveRsity 1.9 (Keenan et al., 2013). Genetic structure was analyzed at multiple spatial scales with STRUCTURE 2.3 (Pritchard Lab, Stanford University; http://pritchardlab.stanford.edu/structure.html) without prior population assignments. A range of assumed population cluster size (K, set sequentially from 1 to 8–12, depending on the analyses) was run 5 times using a burn-in of 500 000 iterations and a run-length of 1 × 10⁶ iterations, using the ΔK criterion (Evanno et al., 2005) to assist selection of “best” Ks. Structure analyses were complemented with factorial correspondence analyses (FCA) implemented in GENETIX.

The genealogic relationships among mtIGS haplotypes were inferred with Network 5 (Fluxus-Engineering; http://www.fluxus-engineering.com/sharenet.htm) using the median-joining algorithm, and their geographic distribution and frequency mapped. Haplotype (H _hap)and nucleotide (Θ_π) diversities were estimated for each population and selected geographical regions using DNASP 6 (Julio Rozas Lab, Universitat de Barcelona; http://www.ub.edu/dnasp/), the former coding indels as single mutations and the latter discarding indels in pair-wise comparisons only.

Ecological niche models (ENMs) were developed with Boosted Regression Trees (BRT; Elith et al., 2008), a high-performance machine learning algorithm that automatically fits complex interactions and nonlinear relationships between predictor variables, while reducing the potential of overfitting by forcing monotonicity responses (Hofner et al., 2011) and by optimal parametrization (Elith et al., 2008). The algorithm fitted georeferenced records of occurrence against biologically meaningful predictor variables for macroalgae species (Fragkopoulou et al., 2022), namely, nutrients (as nitrate and phosphate), salinity, sea ice thickness and maximum winter and summer sea surface temperatures.

Occurrence records of F. distichus were collated from the fine-tuned dataset of marine forests (Assis et al., 2020) which integrates data from different sources (e.g. databases, literature and herbaria) and uses a flagging system to exclude potentially erroneous entries (e.g. records outside known biogeographical ranges). Because absence records were unavailable at the scale of the study, pseudo-absences were randomly generated using a 1:1 ratio with presence records, in locations with no presences recorded (Barbet-Massin et al., 2012) and limited to the biogeographical provinces (Spalding et al., 2007) where the species is known to occur, as well as their neighboring provinces. This allowed restricting pseudo-absences to regions where dispersal can potentially occur, a crucial step of ENM (Martins et al., 2021). The negative effects of spatial autocorrelation and sampling bias (e.g. regional differences in density of records) were further reduced by randomly choosing one record only from those at distances considered significantly autocorrelated, as estimated by a correlogram testing the spatial variability of predictors as a function of geographic distance (Boavida et al., 2016). Predictor variables describing sea surface conditions were obtained from Bio-ORACLE v2.1 (Assis et al., 2018c; Tyberghein et al., 2012) for near present-day conditions and for the Last Glacial Maximum (LGM).

A cross-validation (CV) framework using 10-fold latitudinal bands was implemented to tune the optimal parameter combination of BRT and to estimate the performance of the models and their potential transferability in independent data (Assis et al., 2018a; Krause-Jensen et al., 2020; Fragkopoulou et al., 2021). The “grid search” statistical method interactively tested all possible numbers of trees (range 50–1000, step 50), tree complexity (range 1–6, step 1) and learning rate (0.01, 0.005 and 0.001), by fitting occurrence records with all but one latitudinal band at a time. At each step, performance was estimated with the area under the curve (AUC). The best combination of BRT parameters was identified as that resulting in models with higher AUC in CV (Assis et al., 2018a; Krause-Jensen et al., 2020; Fragkopoulou et al., 2021). Overfitting was further controlled by forcing monotonic responses based on the expected biological effect of predictors (positive or negative) on the probability of occurrence (Hofner et al., 2011). Maximum temperatures and ice fitted negative monotonic responses, whereas nutrients and salinity fitted positive responses (e.g. Assis et al., 2018a; Fragkopoulou et al., 2021; Martins et al., 2021; Krause-Jensen et al., 2020).

ENM were performed using two approaches: (1) species modelling considering the entire amphiboreal range and (2) sub-taxon modelling (i.e. within-taxon niche structure; Assis et al., 2016), using distinct models for the Pacific and Arctic/Atlantic ranges where differences in evolutionary history, genetic composition and diversity were apparent (see results and discussion sections). Habitat suitability maps for the present and the LGM were produced with a BRT model using the optimal parameters found in CV. These maps were further reclassified to presence and absence (binomial responses) using a threshold maximizing specificity (true negative rate) and sensitivity (true positive rate; Jiménez-Valverde and Lobo, 2007).

Niche overlap between the Pacific and Arctic/Atlantic ranges was estimated with the Bayesian framework of Swanson et al. (2015). This defines the niche of each range as the 95% probability region of the multivariate environmental space defined by ENM, and determined overlap as the probability of an individual from the Pacific being found in the niche of Arctic/Atlantic, and vice-versa. Additionally, physiological tolerance limits for the Pacific and Arctic/Atlantic ranges were captured from partial dependence functions (Elith et al., 2008).

A total of 31 collections of Fucus distichus (590 individuals) were sampled from 16 Pacific and 10 Arctic/Atlantic sites ( Table 1 ), covering most of the species range in the northern hemisphere. Collections from Alaska (Pirates, Kayak) and Oregon (Bob Creek) contained a few F. spiralis L. that were easily recognizable genetically and that were excluded from analyses. Microsatellite polymorphism (f_allele > 0.005) ranged between 2 and 13 alleles per locus, averaging 6.14. A total of 38 mtIGS haplotypes were detected among the 470 individuals sequenced ( Table 1 ; Figure 1 ; Supplementary Table S2 ). A subset of 16 haplotypes (≈42%) accounted for almost 90% of the sequences, while another 14 were sampled only once. Several indels of variable size were apparent ( Figure 1 ). Genetic diversity within populations was typically low (A ≤ 3.2; H _E ≤ 0.41, mean n _hap ≈ 2.15), as most populations were monomorphic (or nearly fixed) at several loci ( Supplementary Figure S1 , Supplementary Table S2 ). Nearly all populations exhibited significant heterozygote deficits. Avacha Bay (pop. 23), the only populations showing fixed heterozygosity at locus Fg25 (and hence anomalous positive F _IS values), also exhibited significant heterozygote deficit when the putative duplicated locus was removed from analyses (data not shown).

F. distichus exhibited strong phylogeographic structure. MtIGS haplotypes (Genbank accession numbers PP358930–PP359399) formed a complex network of shallow haplotypic radiations corresponding to at least 10 major genetic/geographic discontinuities ( Figures 1 , 2 ). Eight of these were endemic to the Pacific and 2 were endemic to the Arctic/Atlantic, i.e. they were never found together. Atlantic haplogroups were genealogically close and they were not more differentiated from Pacific haplogroups than different Pacific haplogroups were between them. Haplogroup DO was widespread throughout both sides of the Atlantic and in the Arctic (encompassing samples originally assigned to F. distichus s.s., F. anceps, F. edentatus and F. evanescens), whereas haplogroup LO was restricted to the cold-temperate NW Atlantic (Maine and Newfoundland, Figure 2A ). In Maine (populations 4 and 5), these haplogroups matched two regionally parapatric eco-types – LO corresponding to smaller forms from high-intertidal rock-pools (F. distichus s.s.-like) and DO to larger forms from the low intertidal (F. evanescens-like, Figure 2B , top right). In the NW Pacific and western Aleutians (Adak, pop. 22) only the P, DY (endemic) and Y (shared with Alaska) haplogroups were sampled. All 9 haplotypes detected there were endemic to the region and the majority was restricted to specific populations ( Table 1 ; Supplementary Table S2 ). The remaining 5 haplogroups were restricted to the NE Pacific ( Figure 2A ). The LB, DB and DG haplogroups were present in the Mendocinian and Alaska regions, extending into the eastern Aleutians (Unalaska, pop. 21). In the Mendocinian region each haplogroup dominated successively short stretches of coastline, whereas in SE Alaska they co-occurred more extensively (together with the Y haplogroup) at local scales. Remarkably, haplotypic composition among Alaskan populations varied considerably at very local scales, such as among High- vs Low-shore collections at Pirates ( Figure 2B , bottom centre). Finally, the R and LG haplogroups were endemic to the Montereyan region. In Pigeon Point (pop. 11), these haplogroups matched a priori “large” and “small” morphological assignments, respectively.

Despite the relatively small global allelic diversity, the Structure analyses separated Arctic/Atlantic and Pacific genotypic clusters ( Figure 2B , top) with good support (Δk = 2922; Supplementary Figure S2 ). This high-level subdivision was also recovered by the FCA ( Supplementary Figure S3 ). Within the Pacific there was no support for large-scale regional subdivisions (maximum Δk < 30; Supplementary Figure S2 ), to some extent due to low microsatellite allelic diversity and resolution. Finer-scale regional sub-structuring was evident nevertheless when different regions were analyzed separately ( Figure 2B , bottom). In the NW Pacific and W Aleutians, where most populations exhibited fixed allelic differences in at least one locus ( Supplementary Figure S1 ), recovered genotypic clusters (K = 6, Δk = 254; Supplementary Figure S2 ) matched closely single populations ( Figure 2B , bottom left). The same applied from California to British Columbia, where genotypic clusters (K = 6, Δk = 370; Supplementary Figure S2 ) again matched to a large degree specific populations or contiguous population sets ( Figure 2B , bottom right). Despite local sympatry, “large” and “small” morphotypes from Pigeon Point, characterized by different R and LG haplotypes (see above), also matched distinct genotypic groups. Conversely, in SE Alaska, support for 3-5 genotypic clusters was weak (60 < Δk <90, Supplementary Figure S2 ), and these did not consistently match distinct haplogroups, shore-heights, or regional populations ( Figure 2B , bottom centre). Structure analyses also failed to resolve the two eco-morphotypes/haplogroups from Maine (or any other set of Atlantic populations), even when analyzed alone (Δk < 100, Supplementary Figure S2 ).

Genetic diversity was much higher in the Pacific than in the Atlantic for both microsatellite (A ₍₂₃₄₎ = 6.96 vs 3.86) and haplotypic (n _hap = 29 vs 9, Θ_π × 10⁵ = 857 vs 311) data ( Table 1 ). The Atlantic was also less diverse genetically than any of the Pacific regions considered ( Figures 3A, C ; Supplementary Table S3 ). The region between Point Conception (California) and Vancouver (British Columbia), encompassing the Montereyan and Mendocinian biogeographical regions, was (marginally) the most diverse region. Notwithstanding the occurrence of multiple, some quite divergent (e.g. R phylogroup) genetic groups, intra-population diversities were low. Pairwise differentiation, on the other hand, was often large, highlighting the depth of fine-scale spatial structure. Regional diversity in SE Alaska was not obviously larger than in the NW Pacific and western Aleutians, but consistently contained the most diverse populations in the Pacific ( Figures 3A–C ). Diversity was estimated for high and low shore heights at some sites (Pirates, Point Louise), thus underestimating intra-population (=site) diversity. In any case, small D values indicated relatively weak pair-wise differentiation between these diverse populations ( Figure 3D ). Haplotypic diversity was also large in this region, despite lacking unique haplogroups.

We collated 1283 unique occurrence records for the Pacific and 1015 for the Arctic/Atlantic. The ENMs retrieved high accuracy scores in independent data (i.e. in cross-validation procedure). Specifically, the model developed for the entire amphiboreal range achieved an average AUC of 0.829 ± 0.077, while the models for the Pacific and Arctic/Atlantic achieved an average AUC of 0.829 ± 0.077 and 0.848 ± 0.13, respectively. The posterior distribution of the probabilistic niche overlap metric between the Atlantic and Pacific groups was 96.24%, while the probability between the Pacific and Atlantic groups was 93.28% ( Supplementary Figure S4 ). Despite the large niche overlap between ranges, physiological limiting points inferred from partial dependency functions differed between groups ( Supplementary Table S4 ), with the niche of the Pacific showing greater thermal tolerance to maximum winter and summer temperatures (1.61°C and 1.86°C, respectively; Supplementary Table S4 ; Figure 4 ).

In the Pacific, ENMs over-estimated the distribution of F. distichus beyond known range limits in Point Conception (California) and Peter the Great Bay (Russia’s Far East), and polewards beyond the Bering Strait into the Chukchi and Beaufort Seas ( Figure 5B ). In the Arctic/Atlantic, the occurrence of the species was under-estimated in the Canadian Arctic in all models ( Figure 5B ; Supplementary Figure 5 ). The major mismatch however was in the NE Atlantic, where Ensemble and Pacific ENMs predicted a much more consolidated presence in the British Isles and contiguous French regions of Brittany and Normandy, as well as an isolated pocket in NW Iberia ( Figure 5B ; Supplementary Figure S5B ). The Atlantic ENM over-predicted the occurrence of F. distichus only as far south as SW England (Cornwall; Supplementary Figure S5A ). Hindcast projections supported a massive latitudinal shift in the potential distribution of F. distichus in the Atlantic between the LGM and today ( Figure 5A ). Suitable climatic conditions in the NW Atlantic would have been restricted to ice-free regions off Virginia to New Jersey and along the then-emerged Grand Banks of Newfoundland. In the NE Atlantic, projections indicated a much more southern and compressed distribution, between western Iberia and the periglacial palaeo-shorelines of the Celtic Sea and western British Isles. Globally, overlap between estimated local LGM and observed present-day distributions was quite restricted, as most of the modern range of F. distichus was either glaciated, emerged, or both during the LGM. According to models, latitudinal shifts were much less extensive in the North Pacific. During full glacial conditions, F. distichus could have spread from the Korean peninsula and Japan in the western Pacific, and Baja California (Mexico) in the NE Pacific, as far north as the paleo-shorelines of the Beringia Land Bridge. This massive range was probably interrupted to some extent between the eastern Aleutians and Vancouver Island (British Columbia), where the Cordilleran ice-sheet reached the coast. Globally, and unlike in the Atlantic, F. distichus still occurs throughout most of its immense glacial range. On the other hand, post-glacial poleward expansions, e.g. past the Bering Strait into the Arctic, have been much more modest.

Globally, F. distichus encompasses multiple phylogeographic lineages and microsatellite genotypic clusters within which shallower genetic discontinuities were also pervasive at more local scales. Nested patterns of structuration are common among fucoid seaweeds and other low-dispersal marine organisms (e.g. kelp, direct-developing invertebrates) and typically result from complex, non-equilibrium demographic/genetic processes (e.g. contraction/expansion and isolation/secondary contact dynamics) most commonly associated with Pleistocene climatic shifts (Marko et al., 2010; Neiva et al., 2016; Assis et al., 2018b; Grant and Chenoweth, 2018, 2021; Grant and Bringloe, 2020). More striking is the unprecedented level of genetic diversity uncovered. F. distichus produced the most complex intra-specific mtIGS network recovered so far among fucoid seaweeds (Neiva et al., 2010, 2014, 2017), with haplotypic diversity comparable to that of an entire Atlantic radiation comprising 5+ species, including F. vesiculosus, F. spiralis, F. guiryi s.l (i.e. F. limitaneus and F. macroguiryi; see Almeida et al., 2022), F. virsoides and a mtDNA-introgressed lineage of F. ceranoides ( Supplementary Figure S6 ). Fucus distichus shares its diplontic life-cycle, hermaphroditic mating system, intertidal habitats and dispersal ecology with most of these seaweeds, and was similarly severely impacted by Pleistocene climatic shifts. The exceptional phylogeographic diversity of F. distichus therefore seems to be best explained by its distinctive amphiboreal range, as large stable ranges (particularly in the Pacific, see below) offer more opportunities for allopatric divergence and lineage persistence. The greater phylogeographic diversity of circum-boreal Saccharina latissima s.l. vs NW Pacific endemic S. japonica supports this prediction (Zhang et al., 2015; Neiva et al., 2018), and additional studies may confirm the generality of the pattern.

Phylogeographic and genotypic diversities were extremely skewed towards the Pacific, supporting the view that F. distichus diversified primarily in the Pacific and that its Arctic/Atlantic range was colonized more recently (e.g. Coyer et al., 2011b). Different haplogroup compositions indicate nevertheless that Atlantic colonization pre-dates the last (de)glaciation, which is consistent with patterns of (nuclear) trans-Arctic divergence (Cánovas et al., 2011). It should be noted that Pacific and Arctic/Atlantic haplogroups were never found co-occurring throughout the Arctic basin, including along known dispersal pathways in the Canadian Arctic and eastern Greenland (Bringloe and Saunders, 2018; Neiva et al., 2018). Unlike the Atlantic side, where F. distichus post-glacial expansions reached the Canadian Arctic and Russia’s Novaya Zemlya archipelago (e.g. Lobus and Udalov, 2021, the latter not genetically confirmed), lack of reports suggest a much more limited post-glacial Pacific expansion past the Bering Strait into the adjacent Laptev, East Siberian and Beauford seas (Wulff et al., 2009; Wilce and Dunton, 2014). Limited availability of intertidal hard substrata, low summer temperatures (<5°C) and reduced salinity concentrations in some Arctic seas may help explain this asymmetry (Jueterbock et al., 2016), but the limiting factors (and actual extent) of post-glacial Arctic colonization from the Pacific should be further clarified.

It is interesting to note that most genera of brown macroalgae with amphiboreal representatives, including several large brown (e.g. Chorda, Alaria and Saccharina) and Arctic-endemic red seaweeds, are better represented (i.e. have more species) in the Pacific than in the Atlantic (Lüning, 1990; Lane et al., 2007; Sasaki and Kawai, 2007; Adey et al., 2008; Neiva et al., 2018; Bringloe et al., 2021), implying a predominant track of Pacific to Atlantic invasion. The biogeography of Fucus, however, does not suggest a primary Pacific origin. The sister genus Pelvetiopsis is composed by 4 NE Pacific-endemic species (Neiva et al., 2017), whereas all 7+ extant congeners of F. distichus are Atlantic endemics, including closely related F. serratus (F. spiralis is regarded as a recent introduction to the NE pacific (Coyer et al., 2006). The split and ensuing vicariant diversification of Fucus and Pelvetiopsis has been hypothesized to result from a trans-Arctic crossing by a common ancestor following the first opening of the Bering Strait at the Miocene/Pliocene transition ca. 5.6 Ma ago (Cánovas et al., 2011; Sousa et al., 2019). The secondary colonization of the Pacific by a F. distichus ancestor represents the most parsimonious hypothesis (Lindstrom, 2001; Cánovas et al., 2011) but is not consensual (Coyer et al., 2006; Laughinghouse et al., 2015). Alternative scenarios include a trans-Arctic split between main Fucus clades I (F. distichus/F. serratus) and II (F. vesiculosus et al.), or an earlier Pacific colonization by a F. distichus/F. serratus ancestor, in both cases involving the Atlantic colonization by (and Pacific extinction of) the F. serratus ancestor. While available data cannot answer this question, it is clear that the inferred colonization of the Atlantic represents only the latest among other trans-Arctic migrations in the evolutionary history of Fucus.

Another factor that likely contributed to the greater genetic diversity and differentiation of F. distichus in the North Pacific is the much wider ice-free glacial ranges, and the much greater glacial/interglacial range overlap, when compared with the Arctic/Atlantic. ENMs revealed extensive ice-free habitats suitable for F. distichus in the N Pacific during past major glaciations, from the Korean and Baja California Peninsulas to as far north as the palaeo-shorelines of Beringia. Even in the NE Pacific, where the Cordilleran ice-sheet spread along the coast and potentially created a distributional gap from British Columbia to the eastern Aleutians, suitable conditions existed from Vancouver to Point Conception (and beyond), where the endemic (and seemingly reproductively isolated) lineages R and LG were identified. Globally, and notwithstanding major climatic shifts across the last glacial/interglacial transition, available evidence supports long-term persistence of F. distichus throughout much of the core (and very wide) contemporary Pacific range. By comparison, the core distribution of F. distichus in the Atlantic overlaps extensively with regions that were under the direct influence of the LGM ice-sheets and available for colonization only at the onset of the deglaciation. The existence of a lineage unique to the NW Atlantic is in line with ENMs supporting suitable conditions just south of the Laurentide ice-sheet and along the ice-free palaeo-shorelines of the Grand Banks of Newfoundland (see also Assis et al., 2018a; Bringloe and Saunders, 2018). In the NE Atlantic, ENMs also revealed suitable habitat beyond ice-sheets from north-western British Isles to at least SW Iberia, but this putative southern glacial range is only marginally occupied today. The proximity of the glacial leading edge (just south of ice-sheet limits) and the modern trailing edge (western Ireland) necessarily translates into minimal areas of overlap (i.e. of long-term persistence), where genetic diversity could have accumulated and persisted until today.

Oddly, ENMs reveal suitable present-day climatic conditions throughout much of the British Isles and even Brittany, with Pacific models extending potential distributions as far south as northwest Iberia, i.e. much beyond its core distribution in Norway and its absolute modern range limit in Ireland. Even throughout its southern range in the British Isles, F. distichus populations are exceptionally rare and localized, and exhibit much less morphological and ecological plasticity when compared to other N Atlantic regions (e.g. Powell, 1963). Indeed, the single ecotype present – dwarf F. distichus var. anceps – is restricted to specific and disconnected sites along the western-facing coasts of Scotland and Ireland, subjected to extreme and persistent wave action, where other intertidal fucoids are scarce or absent. The species is ubiquitous throughout equivalent cold-temperate latitudes in the NE Pacific, where it nearly reaches the Southern California Bight warm-temperate eco-region (sensu Spalding et al., 2007). Most Arctic-cold temperate seaweeds with amphiboreal or amphi-Atlantic distributions (including foundational kelps such as Laminaria digitata, Alaria esculenta and Saccharina latissima, intertidal fucoids such as Ascophyllum nodosum and Fucus serratus, and red seaweeds such as Phycodrys rubens and Delesseria sanguinea) are still common in the Brittany-Ireland transition zone, and some even have relic populations in upwelling areas of NW Iberia. Most display unique gene-pools and diversity hotspots in these areas, indicating stable presence across climatic cycles (Wares and Cunningham, 2001; Maggs et al., 2008; Waltari and Hickerson, 2013).

The fragmented distribution of F. distichus in Europe could be explained by ongoing expansion from the NW Atlantic. This hypothesis is unconventional since most post-glacial trans-Atlantic migration has occurred, because of more benign glacial conditions, from East to West, but it explains some of the observations above. The steady expansion of a sheltered form of F. distichus (F. distichus var. evanescens) along the Skagerrak and Kattegat to at least Rostock in the western Baltic Sea (ca. 54° Lat.) since its mid-1890 introduction to Oslofjord (eastern Norway, presumably from northern Europe) is also compatible with this hypothesis, as it confirms suitable thermal regimes and potential for southern range expansions (Coyer et al., 2002; Dietrich and Hendrik, 2017). Atypically high rates of hybridization with F. serratus further suggest that contact between the species in this area is evolutionarily recent (Hoarau et al., 2015). Disentangling among these hypotheses – relic populations left after post-glacial contraction vs ongoing trans-Atlantic (W → E) post-glacial expansion – requires assessing whether genetic variation in the NE Atlantic is to some extent unique or chiefly a subset embedded within NW Atlantic diversity. The NE Atlantic populations, including Irish and Scottish ones, are not particularly remarkable for the studied markers, but more variable markers may be necessary to conclusively answer this question.

Perhaps other factors are playing a role. Whether pre-glacial or more recently, F. distichus may have suffered an important functional bottleneck during the NE Atlantic colonization. Migration corridors along the Canadian Arctic and Greenland require tolerance to extreme Arctic conditions and thus selective filters and gene-surfing along obligatory trans-Arctic pathways could have selected for specialized high-latitude form(s) pre-adapted to colder environments. If so, the regional NE Atlantic “niche” could be a fraction of the global species niche and its “relic” distribution along warmer British Isles (and absence further south) a reflection of it. There, its dependency for extremely exposed habitats could be explained by lower heat/desiccation stress attained in low-shore turbulent environments. The greater thermal tolerance of F. distichus in Pacific vs Atlantic ranges inferred from niche modelling (~1.8°C) is certainly compatible with this hypothesis, but empirical evidence for less thermal tolerance in the NE Atlantic range (or larger thermal tolerance in southern populations of the E Pacific) is missing entirely. Another potential factor to consider is inter-specific competition. F. distichus may be outcompeted by other intertidal fucoids toward temperate latitudes, becoming restricted to extremely exposed habitats (British Isles) or nutrient-enriched harbors (Kattegat and Western Baltic Sea) where other intertidal fucoids (and particularly low-intertidal Fucus serratus, see Jueterbock et al., 2016) are scarce or cannot establish. Intertidal fucoid communities are also represented in the southernmost areas of the NW Atlantic and NE Pacific (Silvetia, Pelvetiopsis), but are comparatively less diverse and saturated than in the temperate NE Atlantic. These contrasting hypotheses differ conceptually – they either assume fundamental or realized niche differences ( Figure 6 ) – but both are rarely considered in marine phylogeographic literature. Comparative physiological/transcriptomic responses and reciprocal transplants, beyond the scope of this study, are likely to bring new insights, and help explain the nature of these apparent range/niche differences between the NE Atlantic and the remaining range.

In the Pacific, and to a lesser extent in the Atlantic (e.g. Canadian Maritimes), the occurrence of endemic haplotypic radiations provide good indication of glacial range dissections and confirm that some lineages never expanded post-glacially much beyond past refugia. As predicted, the Montereyan region in the southernmost NE Pacific was home to two unique haplogroups (LG and R), the latter the most differentiated within the complex. At higher latitudes, but still beyond the reach of glacial ice, populations tended to be dominated by single haplogroups. For instance, between Vancouver Island and Cape Mendocino, distinct NE Pacific-endemic haplogroups dominated relatively short stretches of shoreline. Populations from Japan and Russia were also dominated by single haplogroups, some endemic to (or nearly so) the NW Pacific (P, DY). Contrastingly, in southern Alaska, NE Pacific haplogroups were locally mixed, also with the trans-Pacific Y haplogroup, resulting in the most haplotypic (and genotypic) diverse populations within F. distichus range. Alaskan populations, even when separately analyzed as high- and low- intertidal bands, also exhibited the largest observed heterozygosity. Ice-sheet reconstructions suggest that the Cordilleran ice-sheet reached the southern Alaskan coastline, presumably creating a distributional gap spanning from the eastern Aleutians to Vancouver Island. This region apparently lacks endemic haplogroups, as the ones co-occurring in the area were also found farther south or farther west. Regional diversity was also not atypically high, at least when comparing to the Montereyan/Mendocinian region, and pair-wise differentiation between populations was quite depressed, indicating weaker genetic structure and lineage sorting when compared to other regions. Finally, local populations were often composed of parapatric bands reminiscent of distinct ecotypes (even if not consistently recovered as separate entities, see next section).

Taken together, these data suggest that the Alaskan Gulf represents a melting pot where unique and more structured southern East and West Pacific sources met. Most Y and DG haplotypes samples were, however, private to this area. It remains a possibility that they are shared with contiguous, unsampled areas beyond the direct influence of the coastal ice-sheet, but if not, they could signal glacial persistence at these higher latitudes as in the kelps Alaria, Hedophyllum nigripes and Saccharina latissima (Grant and Bringloe, 2020; Grant et al., 2020; Grant and Chenoweth, 2021). The mixed evidence is compatible with a more complex scenario were F. distichus survived glacial advances in SE Alaska in a few small pockets of unglaciated shorelines, before invasion from southern and western ranges following ice-retreat – i.e. the area representing both a relic northern refugia and a post-glacial melting pot.

The Aleutian archipelago was poorly sampled but was also expected to bear signatures of recent contact between NW and NE Pacific regions. Central Aleutian Adak showed a more westerner character. This possibly reflects more benign LGM conditions, which presumably allowed survival (or at least earlier colonization) of western vs eastern Aleutians, as the latter was under the influence of the Cordilleran ice-sheet. Additional collections (analysed at a later stage and not included in this study) reveal that eastern Aleutian Islands indeed are genetically closer to the contiguous Alaskan region.

In line with previous studies, molecular data failed to support traditionally recognized morpho-taxa among Atlantic F. distichus, nor, for the most part, subdividing Pacific “Fucus gardneri”. Much of the morphological variation between ecotypes has been shown to be continuous and explained by environmental factors (Powell, 1957; Rice et al., 1985; Sideman and Mathieson, 1985; Ruuskanen and Nappu, 2005), although a few common garden experiments suggest that at least some variation is heritable and not environmentally-induced (Sideman and Mathieson, 1985). In any case, eco-morphological or genetic entities were largely incongruent – some phylogroups included multiple ecotypes, and vice-versa, similar ecotypes belonged to distinct genetic lineages. In other words, at large spatial scales, F. distichus is best described as a single ecologically plastic, morphologically polymorphic and genetically structured species whose variation (genetic, morphologic, etc.) is not sufficiently discontinuous or consistent as to justify recognizing separate genetic or taxonomic entities.

Still, it remains challenging to fully capture and define the biodiversity contained within F. distichus. In some regions clearly distinct ecotypes matching distinct gene-pools co-existed at very local scales as stable sympatric/parapatric entities. This pattern seems to imply (local) isolation mechanisms and is reminiscent of incipient species (Billard et al., 2005; Neiva et al., 2018). Fucus are externally fertilizing broadcast spawners and show broad gametic compatibility, as revealed by natural hybridization/introgression between many species pairs (Coyer et al., 2002; Neiva et al., 2010). Remarkably, hybridization rates between F. distichus and F. serratus are reduced in naturally sympatric ranges when compared to recently established (19^th century) contact zones, implying reinforcement, i.e. the strengthening of pre-zygotic isolating barriers in response to selection against less-fit hybrids (Hoarau et al., 2015). Local build-up of isolation among younger and less differentiated conspecific populations has also been documented in the context of range-shifts, ecological speciation and ring-species (Hopkins, 2013). Many F. distichus haplogroups evolved in allopatry and in some cases may be partially pre-adapted, as suggested in parapatric contexts, to different environmental conditions. Upon secondary contact, their local and regional integrity may be explained by limited gene-flow resulting from high selfing rates, spatial segregation and/or temporal reproductive asynchrony, potentially reinforced by reduced hybrid fitness and/or strong selection for contrasting environmental conditions (Sideman and Mathieson, 1983; Coyer et al., 2002; Billard et al., 2010; Zardi et al., 2011; Hoarau et al., 2015; Monteiro et al., 2016). In Maine, for instance, a high-intertidal tide-pool ecotype (LO haplogroup) reproduces during low-tide when it is isolated from larger, lower-shore (DO haplotype) conspecifics (Pearson and Brawley, 1996), but low microsatellite resolution prevents at present confirming their degree of reproductive isolation. Dwarf and normal-sized intertidal ecotypes from Maine have also been shown to exhibit reproductive asynchrony (Sideman and Mathieson, 1983). In Alaska, high- and low-intertidal bands often showed different haplotypic compositions, but the pattern tended to be site-specific and only weakly correlated with respective nuclear backgrounds, requiring further examination. This area was also the one where populations exhibited the largest observed heterozygosity, implying an important incidence of outcrossing. The most convincing example is perhaps that of Pigeon Point, where haplogroups R (wide-frond type) and LG (thin-frond type) co-occurred sympatrically. Despite overlapping vertical ranges and reproductive period, only a single genetic F1 hybrid was detected (assuming it is not a laboratory contamination during DNA extraction), suggesting nearly complete reproductive isolation. Ongoing analyses of nearby populations support this pattern at a regional scale. This R haplogroup is also the most divergent within the complex, so there is convincing additional ground for suspecting a separate, incipient species.

From a taxonomic perspective, this “paradox” – the existence of putatively isolated genetic entities that are relatively easy to recognize and discriminate at local scales but much harder to define (in terms of distinctive and constant morphological, ecological or genetic characters) at larger spatial scale - is problematic. Dealing with these entities may be even more complicated if reproductive isolation, as seems to be the case in some regions, is contingent and varies geographically depending on the specific genetic entities involved and/or the evolutionary and environmental contexts underlying their regional assembly. Similar studies may show that similar patterns apply to other highly-structured plastic species spanning large geographical areas regularly heavily impacted by glaciations.

Widespread amphiboreal marine species such as F. distichus are ideal models to investigate how climatic shifts, dispersal barriers, environmental filters and distinct biotic assemblages influence intra-specific genetic structure and niche differentiation across species ranges. Despite the unprecedented effort, a range of open questions remain that could be further addressed in future studies. First, new markers and approaches should be employed to estimate the timing of trans-Arctic migrations and secondary colonization of the Atlantic, and particularly of the NE Atlantic rear-edge, that are difficult to establish with microsatellite and mitochondrial data alone. The same applies with respect to spatially co-existing ecotypes, for which genomic data can provide much more power to clarify the actual degree of divergence and isolation and best taxonomic status, including confirming the putative cryptic species in the eastern Pacific rear-edge. As mentioned above, experimental work also seems necessary to clarify the nature of the apparent and very peculiar niche unfilling in the NE Atlantic. Finally, it is unquestionable that ongoing climatic change will translate into wider areas and windows of ice-free conditions in the near future. Deeper expansions into the high-Arctic seem inevitable and might even benefit from increased shipping along Canadian and Russian Arctic shores (Jueterbock et al., 2016). Vicariant Pacific and Atlantic lineages may eventually meet, creating an excellent biogeographic opportunity to directly investigate the fate of divergent lineages upon secondary contact (e.g. Neiva et al., 2018). These opportunities of research will not only improve our understanding of the complex biogeographical history, trans-Arctic divergence, and species cohesion in Fucus distichus, but of amphiboreal communities in general.