A literature review located the historical southern endpoint of F. vesiculosus in Southern Morocco at Khnifiss Lagoon (27°59'40.82"N-12°16'33.98"W; Figure 1). Its presence at this location was last reported in 1985-1986 and it was also confirmed by the Institut Scientifique, Département de Botanique et d'Ecologie Végétale in Rabat (personal communication by Dr. Ibn Tattou Mohamed; see Additional file 1).

The southern endpoint is now positioned on the Portuguese west coast in the Tejo estuary (38°45'38.99"N, 8°56'28.43"W), representing an 11° latitude shift (1,250 km range contraction; Figure 2b). In addition, an isolated population was located in Cadiz-Puerto Real (36°31'27.53"N, 6°10'52.61"W) approximately in the middle of this range contraction and unattached vegetative fragments of adult F. vesiculosus were still observed entangled in Spartina grass in Ria Formosa (37°00'36"N, 7°59'31"W), although the previously permanent and attached populations have totally disappeared from this coastal lagoon.

Historically, F. vesiculosus populations were reported along the whole Northern Iberian coastline, with the latest presences at north eastern locations dated from 1951 to 1988. Currently, a limited retreat of 70 km within the Bay of Biscay has occurred and the present endpoint is at Colindres, eastern Cantabria (43°23'36.63"N, 3°27'50.15"W).

Temperature isotherms shifted polewards to higher latitudes (Figure 2c). This gradient is interrupted by areas of strong and persistent upwelling. In particular, in Morocco, the constant upwelling off Cape Ghir (centered at 31°N) generates a clear cold thermal intermission. Other upwelling cells are present to the south (Western Sahara), but fall outside the distributional range of F. vesiculosus. In contrast to North African shores, upwelling waters off north-west Iberian shores contribute to the SST gradient.

Over the last 30 years, a significant increase of coastal SST (°C/decade; Figure 2a) was observed throughout most of the entire Iberian and North African Atlantic coast. Coastal warming rates were not significant at relatively restricted areas within upwelling cells off Morocco (between 31°22' and 32°22') and Western Sahara (between 24°22' and 25°52'). In non-upwelling areas, coastal SST warming ranged from 0.092 to 0.297°C/decade, with an average of 0.214°C.

Previously sampled F. vesiculosus collected in 2001 and 2002 from 12 locations, including four extinct populations (Figure 1), were genotyped at five microsatellite loci. The whole dataset (668 individuals genotyped) was highly polymorphic and the total numbers of alleles per locus were: 17 (locus 20; Additional file 2), 8 (locus 58; Additional file 3), 11 (locus 94; Additional file 4), 14 (locus 78; Additional file 5) and 9 (locus 38; Additional file 6).

Allelic richness (Â) varied between 6.28 and 3.08 with lowest values for ML, RF, TV and LX, while unbiased heterozygosity (H_E) ranged between 0.381 and 0.587 with highest gene diversity values for RE, VN, RL and RM (location codes as reported in Figure 1). Observed heterozygosity (H_O) was significantly different from H_E in nine locations where heterozygote deficiency was detected, resulting in significant positive values of the inbreeding coefficient F_IS (Table 1). The high values of F_IS obtained were not locus dependent, indicating the markers used did not display technical issues such as null alleles. Unique alleles were detected for 8 of the 12 populations and varied between one and three.

Strong differentiation among populations was found as estimated by pairwise genetic differentiations F_ST and D_EST values (10,000 permutations, P < 0.01 for all comparisons, see Additional file 7; all confidence intervals do not overlap with zero, see Additional file 8). Pairwise differentiation was lower within two subgroups (hereafter 'northern' group: OV, RE, LG, VN, RL, ML; and 'southern' group: MG, TJ, RM, RF, TV, LX; location codes as reported in Figure 1). However, D_EST estimates were particularly high between RE and all the other locations whereas LX and TJ had the lowest divergence.

On the basis of the dimensions scores of the Correspondence Analyses (CA), populations were grouped into two major clusters corresponding to the northern and the southern subgroups (Figure 3a). The first axis explains mostly the variance between these two clades (34.86%), while the second and third axes mainly illustrate the differentiation within the two subgroups (13.48% and 11.77%, respectively). Overall, the northern group had a higher inter-population diversity, which is illustrated by the more scattered distribution of the points in the plot. The CA is consistent with the tree topology based on allele frequencies using Cavalli Sforza and Edwards' chord distance (Figure 3b), where two main branches were resolved. The first of these two subgroups comprised six populations (northern group). A second subgroup (southern group) consisted of six populations, two of which are still present at the locations while the rest are now extinct. This major genetic subdivision was supported by the results obtained with STRUCTURE, which revealed the most significant increase at ΔK = 2 clusters (Figure 4a), thereafter ΔK remained unchanged. We concluded that K = 2 is the most likely number of genetic clusters (Figure 4b).

The northern group had many unique alleles, but the southern cluster including the now extinct populations also displayed unique alleles. This cluster suffered diversity loss with the recent extinctions, but still retains unique diversity (three unique alleles detected in five loci assessed, at MG, TJ; Table 2).

Systematic shifts in geographical ranges of species and increased extinction risk has emerged as one of the most pervasive biological responses to recent anthropogenic climate warming. The general trend shows range expansion towards higher latitudes; however, redistribution dynamics at southern range margins are less clear (for example, 5). While on land, southern or 'trailing' range boundaries of terrestrial ectotherms may not shift consistently towards the poles with climate warming, marine species ranges match more closely their limits of thermal tolerance. Thus, range shifts will show a more predictable southern range contraction in the marine realm 4.

Extinction of southern edge populations along the Iberian Peninsula has been reported for other marine cold temperate species 17. Another study conducted along this stretch of coast highlights the complexity in making generalizations about range shifts for intertidal species 37. The extensive historical records reviewed in our study indicate that until the late 1980s abundant populations of F. vesiculosus were still present along the entire Moroccan and Iberian Atlantic coasts, while in the early 2000s extensive populations of this alga still persisted in southern Portugal. When these records are compared with our more recent surveys (2009 to 2011), a major contraction of the southern range margin of this keystone intertidal species (that is, from southern Morocco to central Portugal, latitudinal reduction of more than 1,200 km) is revealed. Taken together with past studies, there is a strong indication that the northward distributional shift along the Portuguese coast occurred very recently, probably within the last decade. This range contraction is an order of magnitude larger than the range shifts recorded for other species over the last 50 years along the Northeast Atlantic (for example, 3137).

Trends in SST indicated that the large scale disappearance of southern populations of F. vesiculosus is correlated with warming of North African and southwest Iberian coastal waters, with steady year-round increases averaging >0.2°C/decade over the last three decades. Seasonal thermal extremes might also play a crucial role in defining the distribution of species through sudden population declines or die-offs 39. However, the key pressure exerted by increasing average SST may result from sub-lethal rather than catastrophic effects. In the North-Eastern Atlantic, the distributional limits of most warm and cold-water algal species seem to be set by sub-lethal effects of water temperature, through reduced reproduction and growth, rather than by lethal effects 3440. Several other factors indirectly related or unrelated to changes in temperature may be responsible for the observed distributional shifts. These include biotic interactions (for example, grazing and competition 41) and near-shore abiotic or anthropogenic effects (for example, pollution, coastal erosion, wave 3842). Moreover, seaweed survival, growth, and reproduction are known to be largely affected by latitudinal variations in nutrient supply 4344. However, mean or maximum water temperatures appear to be the best overall predictors of mortality events in the intertidal zone and the main determinants of large scale range contractions (for example, 36454647) as opposed to maximum air temperatures. Intertidal organisms are essentially at the same temperature as water during high tide, while thermal conditions experienced during low tide are more likely to be affected by individual physiological 4849 and behavioral 50 capacities and by the heterogeneity of intertidal microhabitats 5152.

The abundance of species is generally expected to peak at the center of their distributional range and decline with adverse conditions towards the range limits (for example, 53). However, several empirical studies have challenged the 'abundant center' view of species distribution 54, such as in the intertidal zone where localized variations in the thermal environment due to climatic and tidal interactions disrupt monotonic latitudinal temperature clines 545556. In particular, delivery of cold waters by upwelling is not directly correlated with atmospheric processes, providing an opportunity for regional de-coupling of global warming and the creation of refugial areas where intertidal communities can escape the effects of projected climatic changes 5758. Strong upwelling events off the northwest Iberian coastline reinforce latitudinal temperature gradients 32, while local upwelling further south (for example, the Cape Ghir region in Morocco, 31°N) disrupts latitudinal thermal clines. The latter constitutes one of the most persistent (year-round) upwelling cells along the northwest African coast 59. The SST data clearly indicate that the upwelling region around Cape Ghir has not experienced significant warming, in contrast to adjacent non-upwelling areas (but see below). As a result, and in contrast to the Portuguese coastline, this upwelling cell disrupts the regional thermal latitudinal cline. Despite suggestions that upwelling could provide refugia from the thermal stress of ocean warming 5758, our results show that F. vesiculosus populations affected by cold upwelled waters have also undergone extinction. In Cape Ghir, warming rates show variation at monthly timescales, with a significant SST increase of about 0.28°C/decade in March (http://coastalwarming.com 60) and maximum SST peaking occasionally over 24°C 61. Although monthly SST anomalies and high temperature events are less frequent and intense than in contiguous non-upwelling areas, they might be strong enough to cause extinction. Furthermore, effects of climate variability can be aggravated if isolation limits gene flow between populations, increasing genetic loss to drift with consequent loss of evolvability to cope with future environmental change 9. While northern populations are well connected by continuous patches of F. vesiculosus, populations at the southern edge were fragmented and spatially isolated prior to extinction, potentially limiting population size and increasing susceptibility to drastic demographic and genetic variation.

Interestingly, patches of F. vesiculosus persist within the Bay of Cadiz, along the retreating front. This represents a restricted (extending less than 1 km along the coast) population confined to very sheltered areas within the bay (personal observation). The forecasted warming trend suggests that this small and isolated relict population faces a high risk of extinction, and this prediction is further supported by the extinction of nearby (less than 10 km distant) populations (for example, La Caleta and Playa de la Victoria).

Along eastern Atlantic shores, F. vesiculosus extends its northern distributional limit in Norway, also inhabiting the northernmost and brackish regions of the Baltic Sea and in the White Sea into brackish tidal marshlands. Within the range of F. vesiculosus, there is a well-documented gap along the SW French coast (French Basque Country and Gascony) due to lack of rocky substratum 62. In contrast to the southern clade, our results indicate that the present distributional range of the northern lineage in Iberia has been affected minimally, with a contraction of approximately 70 km at the southern limit of this gap (Spanish Basque Country). Notably, F. vesiculosus individuals from other regions in Europe are genetically closer to the northern Iberian clade identified in our study by the same markers 63. A graphical model based on IPCC (Intergovernmental Panel on Climate Change; 64) scenarios predicted the range decline observed in our study 65. This model extends its forecast until 2050 predicting that, by 2025, cold water species, including F. vesiculosus, will disappear completely from the Bay of Biscay.

In Northern Spain, recent distributional declines have also been reported for Fucus serratus 66, Himanthalia elongata 67 and three kelp species 68, most probably caused by significant increases in SST 69 and unusually intense warm inflow of seawater during autumn-winter in 2006 70. Interestingly, despite these warming trends, present SST at the eastern limit of F. vesiculosus in Northern Spain is lower than that at the southern limit in Portugal. This suggests that the two genetic lineages have different thermal tolerances or that distinct factors are setting their distribution limit. Indeed, different genetic lineages or populations inhabiting different areas within a species' distributional range can show diverse resistance or resilience to abiotic stresses (for example, 4971). For example, along Atlantic European shores, southern edge populations of F. serratus are less resilient to desiccation and heat shock than central populations 49. Future experimental approaches could allow explicit tests for signs of local adaptation between southern and northern lineages of F. vesiculosus.

Reduction in genetic variation gives independent evidence of the severity of population and habitat contraction. Estimating losses in genetic diversity associated with habitat contraction and environmental change is one of the key challenges in biodiversity research. The intensity of genetic loss largely depends on the geographical distribution of genetic diversity within the species, including cryptic diversification. Our analyses revealed two clusters, coinciding with more northern and southern locations, distinguished by several analyses based on allele frequencies, and containing unique alleles within each group, in the five microsatellite loci here used. The differentiation between northern and southern lineages in F. vesiculosus was clearly established using a large number of distinct genetic markers including 13 protein coding genes 72 and 35 SNP markers 73. The separation between these lineages of F. vesiculosus is older than the divergence between other more recently evolved species within the genus Fucus 71, highlighting their conservation value. The same multigene assays 7071 could not be performed on the old DNA left from the extinct populations here studied, thus our genetic analyses of the extinct populations are only based on five microsatellites, yet these were sufficient to reveal unique alleles.

The divergence between northern and southern lineages does not appear to be absolute, since a few individuals appear to occur in mixed populations (STRUCTURE analysis; for example, RE, OV, RM). Additionally, individuals from RE, not only diverge from southern populations, but emerge as the most distinct within the northern clade (CA analysis) with no apparent geographical explanation. This appears geographically unexpected, but patterns of genetic structure and differentiation are not necessarily maintained by prominent ecological/oceanographic barriers to dispersal 74. For example, in Fucus ceranoides, distant (but relatively similar) populations do not necessarily exchange more migrants than closer (but relatively distinct) populations 18; historical patterns of extinction and colonization may play a more important role than ongoing gene flow in determining the extent of genetic divergence between extant populations.

Presently, the distributional range of the southern clade is restricted within approximately two degrees of latitude (17% of the historical range) and southern intra-specific genetic diversity is represented by few extant populations under threat from climatic warming. Given that predictions of future climate indicate a further rise in SST 75, there is a real risk of extinction of the southern phylogeographic distinctiveness owing to climatic shifts. In contrast, for long-term conservation purposes, it is reasonable to assume that, despite the relatively small range decline, the diversity of the northern clade is not under immediate risk of reduction because populations of this same genetic group are still present throughout Northern Europe.

Past and present distributional ranges of the southern lineage are characterized by markedly higher SST than those of the northern lineage; it can be hypothesized that the existence of this lineage is explained by local adaptation or phenotypic plasticity to withstand warmer waters. Our multi-locus genetic data represent neutral genetic change and do not necessarily reflect evolutionary responses to selection. However, divergent lineages often arise in areas of speciation (for example, 7677). Indeed, recent data indicate that the hermaphroditic species F. virsoides, F. spiralis and the recently-described F. guiryi 78 are derived from a dioecious ancestor sister to the southern clade of F. vesiculosus 72. The potential loss of the cryptic southern F. vesiculosus lineage could end ongoing diversification (or speciation) processes and compromise the adaptive potential of the species as a whole in the face of future global warming.

In this study, genetic diversity of F. vesiculosus, as expressed by both gene diversity and allelic richness, was higher in northern Iberian populations than in the southern lineage. Higher genetic diversity in the northern Iberian region can be a signature of higher temporal stability of large populations while the opposite, in the southern region, could be the result of lower population sizes and/or temporal variability (that is, bottlenecks, extinctions/recolonizations) of the trailing edge populations. Higher diversity could also result from secondary contact between distinct lineages and/or from hybridization and introgression. While F. vesiculosus co-occurs with other fucoids throughout the majority of its range, including Northern Iberia, Southern Iberian individuals occur in allopatry. Therefore, hybridization events are likely only within and to the north of the contact zone between the sympatric and allopatric ranges 63. Given the apparently small population sizes (small patches) of the populations of the southern group, which are, thus, prone to drift and bottlenecks, and their geographical isolation (separated by tens to hundreds of km), it is perhaps surprising that genetic diversity, although lower, still remains so close to northern levels. Regardless of the cause behind the current within-clade diversity, their genetic distinctiveness implies that the total diversity of F. vesiculosus would be significantly reduced should these last remaining populations also become extinct.

Data on the past distribution and southern endpoint of F. vesiculosus along Atlantic North African and Iberian coastlines were gathered from a comprehensive literature review and information from herbarium collections. Data for the last decade also included personal records and associated sample collections made by our team in multiple surveys of the Iberian and Moroccan coastline for genetic sampling of this species for previous studies 63728283 and of species that co-occur in the same habitats 178485. Published literature was screened up until May 2012 in searches using Google Scholar and ISI Web of Knowledge and based on combinations of the following keywords: Fucus, vesiculosus, Morocco, Iberia, Iberian Peninsula, Portugal, Spain, France, Atlantic, Western Sahara, distribution, herbarium, Fucales, Fucophyceae. In addition, literature reported in AlgaeBase (http://www.algaebase.org) was screened and searched up until 16 May 2012. Sampling artifacts may bias range shift studies 86 because (a) the spatial and temporal scale of surveys are limited, and (b) taxonomical information and ecology of species are poorly known. We expect these limitations to have minimal effects on our data because (a) the entire stretch of coast examined in the present study has been the object of extensive research providing a fundamental dataset of the spatial and temporal extent of the species' geographic range; and (b) F. vesiculosus is a scientifically popular species with well-known biology. It has the clear distinctive trait of possessing air bladders on the thallus, diagnostic for this species, and presents separate sexes, in contrast with the only other congeneric species occurring southwards from North Portugal. In addition, museum and literature records did not show any discrepancy. Therefore, we are confident about the quality of the points of reference used in this study.

Current distribution and endpoints were obtained based on an extensive field survey thoroughly covering the entire past distribution of F. vesiculosus. The survey was carried out during low spring tides between 2009 and 2011 on rocky intertidal shores along the Atlantic North African and Iberian coastlines. Eighty-four moderately wave exposed and estuarine locations were sampled from Hondarribia (Spain, 43°21'N 1°47'W) to Dakhla (Western Sahara, 23°41'N 15°55'W), covering a latitudinal extension of more than 2,150 km. Locations of the survey were chosen based on previous reports of the presence of F. vesiculosus (herbarium, DNA samples and literature data). Five additional locations were only surveyed in 2009 to 2011, of which two were in Western Sahara (Boujdour and Dakhla), south of the historical southern limit.

Most of the locations were visited at least twice, covering winter and summer months, and each location included two sites from 500 m to 1 km apart. At each site, two observers performed searches lasting approximately 60 minutes across all microhabitats present.

Monthly averaged SST data between January 2009 and December 2011 with a 4 km resolution were retrieved from the Moderate Resolution Imaging Spectroradiometer-Aqua (MODIS-Aqua) dataset available from the National Aeronautics and Space Administration (NASA) Goddard Earth Sciences (GES) Data and Information Services Center (DISC). Visualization was performed using Giovanni, a web-based application developed by the GES DISC 61.

Coastal SST warming data over the last three decades (°C/decade from January 1982 to December 2011) were obtained from the Worldwide Coastal Warming Assessment project website at http://www.coastalwarming.com/index.html 60. This project uses Advanced Very High Resolution Radiometer (AVHRR), NOAA Optimum Interpolation ¼ Degree Daily Sea Surface Temperature Analysis data (Reynolds OI V2 SST data) acquired from the NOAA-National Climatic Data Center (http://www.ncdc.noaa.gov).

Note that, although warming rates at coastal and ocean locations often depend on the selected spatial and temporal scales, our SST trends are comparable with those reported in other studies that focused on different periods and spatial resolutions 31598788.

Previously sampled F. vesiculosus collected in 2001 and 2002 from 12 locations, including four extinct populations, were used for genetic analyses (Figure 1). DNA extraction used the DNeasy™ Plant Mini Kit (Qiagen, Hilden, Germany) and five polymorphic microsatellite loci, L20, L58, L38, L94, L78 89 were amplified and genotyped following the methods described in Perrin et al. 83.

Observed (H_O) and expected (H_E) heterozygosities were estimated, and deviations from Hardy-Weinberg equilibrium were tested for significance with 10,000 permutations, using GENETIX 4.05 software (Belkhir et al. 1996-2004), applying q-value correction for multiple tests 90. Allelic richness (Â) was estimated for each population separately, standardized to 20 individuals to account for differences in sample size, using the StandArich package and R 2.10.1 software. Plot of allelic frequencies for each location and marker was performed with the above package.

To depict the global genetic variation among the samples, a CA based on the matrix of individual genotypes identified with five nuclear markers was performed using the factorial correspondence analysis (FCA) procedure on populations implemented in GENETIX 4.05 91 and graphically visualized with IBM^® SPSS^® 20 software.

F_ST and confidence intervals were estimated between pairs of populations with the estimator θ 92, and computed using the R package 'DiveRsity' 93. Additionally, pairwise population differentiation was calculated as estimates of Jost's D_EST 94 using SPADE software 95. Significance was tested using 10,000 random permutations of the individuals between samples with a threshold adjusted using q-value correction for multiple comparisons.

The Cavalli-Sforza and Edwards' chord distance 96 was computed (with GENDIST), as this measure has been shown to generate higher probabilities of obtaining the correct tree topology 97. Neighbor-joining was used to assemble the tree in NEIGHBOR with bootstrap re-sampling (10,000 replications) executed using SEQBOOT and CONSENSE. All programs are part of the software package PHYLIP 3.69 98. The output was visualized in FIGTREE 1.3.1. 99 and edited in Adobe^® Illustrator^® CS4 14.0.0 (Adobe System Inc.).

STRUCTURE 2.3.3. software 100 was used to estimate the most probable number of population clusters (K). The analysis was run without prior information on populations, assuming correlated allele frequencies and admixture. The number of possible Ks assessed was 1 to 13 (maximum number of populations plus one) and 20 independent runs with 100,000 Markov Chain Monte Carlo (MCMC) iterations and 10,000 burn-in were performed at each K to calculate ΔK as in Evanno et al. 101. CLUMPP software 102 was subsequently used to find the optimal alignment of the 20 replicate cluster analyses of the same K. The mean membership matrix across replicates was plotted with the program DISTRUCT 103.

Finally, to compare the level of diversity within clusters, estimates of Â, H_O and H_E were calculated for the different clusters determined by STRUCTURE analyses by pooling all populations from each cluster.