The Influence of Pleistocene Climatic Changes and Ocean Currents on the Phylogeography of the Southern African Barnacle, Tetraclita serrata (Thoracica; Cirripedia)

The evolutionary effects of glacial periods are poorly understood for Southern Hemisphere marine intertidal species, particularly obligatory sessile organisms. We examined this by assessing the phylogeographic patterns of the southern African volcano barnacle, Tetraclita serrata, a dominant species on rocky intertidal shores. Restricted gene flow in some geographical areas was hypothesized based on oceanic circulation patterns and known biogeographic regions. Barnacle population genetic structure was investigated using the mitochondrial cytochrome oxidase subunit 1 (COI) region for 410 individuals sampled from 20 localities spanning the South African coast. The mtDNA data were augmented by generating nuclear internal transcribed spacer 1 (ITS1) sequences from a subset of samples. Phylogenetic and population genetic analyses of mitochondrial DNA data reveal two distinct clades with mostly sympatric distributions, whereas nuclear analyses reveal only a single lineage. Shallow, but significant structure (0.0041–0.0065, P<0.01) was detected for the mtDNA data set, with the south-west African region identified as harbouring the highest levels of genetic diversity. Gene flow analyses on the mtDNA data show that individuals sampled in south-western localities experience gene flow primarily in the direction of the Benguela Current, while south and eastern localities experience bi-directional gene flow, suggesting an influence of both the inshore currents and the offshore Agulhas Current in the larval distribution of T. serrata. The mtDNA haplotype network, Bayesian Skyline Plots, mismatch distributions and time since expansion indicate that T. serrata population numbers were not severely affected by the Last Glacial Maximum (LGM), unlike other southern African marine species. The processes resulting in the two morphologically cryptic mtDNA lineages may be the result of a recent historical allopatric event followed by secondary contact or could reflect selective pressures due to differing environmental conditions.


Introduction
Past climatic changes in the form of glaciations, as well as associated changes in temperature and precipitation, have been suggested as drivers of demographic change in both terrestrial and marine species. Numerous marine species have been shown to be affected by sea level fluctuations that occurred during the last 100 000 years and especially around the Last Glacial Maximum (LGM, 26 000 to 19 000 years ago) [1], [2]. Documented changes have included population bottlenecks and alterations of gene flow, as well as extinctions of local populations [3], [4], [5], [6], [7], [8]. In contrast, refugia provided a means of persistence for some populations, especially in species in the northern hemisphere [6], [9], [10], [11], [12], [13], [14]. The effects of sea level changes become less obvious in areas that had no ice cover during glaciation periods, such as southern Africa. Although the documented drop of ,130 m in sea level caused large areas of the Agulhas Bank to emerge with the formation of the southern coastal plain [2], [15], the effects on the topology of the southern African coastline are not well known. Further, no evidence for vicariance events, such as those shown along the south-western coast of Australia as a result of the land bridge across the Bassian Isthmus, [16], [17], [18], have been documented [19]. However, a number of southern African marine species show some degree of population genetic structuring along their geographic ranges, [20], [21], [22], [23], [24], yet barriers vary between species and there are few congruent patterns [19], [25], even in closely related species [26]. The variability in these patterns makes the explanation of the processes involved in generating the population genetic structuring particularly complex [26].
Such variance in patterns also applies to the study of demographic change; some marine species have been shown to maintain demographically stable populations throughout the LGM, without significant changes in population sizes, while others show population bottlenecks with post-LGM recolonization; both scenarios seem equally plausible [8] and can be found even in closely related species [26]. To better understand the evolutionary processes that drive changes in population sizes and distributions, it is important to carefully consider the timing of population expansions and contractions [27]. Further, it has been hypothesized that species with different life-histories show different patterns in response to climatic oscillations [8], [18]. For example, sessile organisms such as suspension feeders may be more likely to persist in cooler climates due to the greater physiological costs and potential harm of mobile lifestyles under harsh conditions [8]. In addition, it has been suggested that the majority of species showing results consistent with LGM persistence are found in the lowest region of the intertidal or shallow subtidal and possess planktonic larvae [8], [18].
In an attempt to contribute to the growing body of literature describing genetic patterns in southern African marine species and to determine the effect of the LGM on an intertidal organism, we chose the volcano barnacle, Tetraclita serrata, one of the dominant invertebrates of the mid intertidal zone with a wide distribution ranging from KwaZulu-Natal (east coast of South Africa), extending into Namibia (west coast of southern Africa; Fig. 1). Although this range spans at least four of the five biogeographic provinces described for the region, the species is more commonly found on the south and east coasts of South Africa ( Fig. 1) [45]. A systematic study, incorporating both molecular and morphological data, suggests that T. serrata consists of two sympatrically distributed populations in South Africa [46]. Some conclusions were made on the demography of the species based on molecular diversity indices, but in the absence of molecular dating, the effect of the LGM could not be established. Furthermore, although hypotheses on gene flow were presented, no gene flow analyses were conducted.
The present study provides a more in-depth phylogeographic consideration, based on a considerably larger data set, to elucidate the intraspecific variation of T. serrata across sampling localities along the South African coastline. As barnacles are obligate sessile organisms with planktotrophic larvae, it is hypothesized that the longevity and dispersal capabilities of juvenile stages and the influence of the associated current systems may dictate the dispersal and genetic connectivity of the species. Given their position in the mid-intertidal, it is also postulated that T. serrata was not severely affected by the LGM (see [8], [18]). To test these hypotheses, the cytochrome oxidase subunit 1 (mtDNA COI) and the nuclear internal transcribed spacer 1 (nDNA ITS1) were used to document the phylogeographic structure of T. serrata; gene flow patterns and demographic changes of the species were investigated using the mtDNA data.

Study area
The oceanographic patterns of southern Africa are complex and dynamic with seasonal effects on offshore, as well as coastal marine species [47]. The area is governed by two major boundary currents; the cold-water Benguela Current, which flows northwards along the west coast of the continent, and the warm-water Agulhas Current, which moves southwards along the east coast of South Africa (Fig. 1). The diverse oceanographic and biological regimes along the South African coastline have been recognized to fall into five distinct biogeographic provinces, each characterized by a unique suite of potentially physiologically-adapted plant and animal assemblages [48].

Collection of material
A total of 410 adult individuals of T. serrata were randomly collected across 20 intertidal sites ( Fig. 1; Table 1). Sampling localities span most of the distribution and have been partitioned according to coastal region (west, south-west, south and east; Table 1), which broadly correspond to four of the five coastal bioregions as defined by Griffiths et al. [48]. All samples were collected in the form of whole individuals and stored in 100% ethanol until DNA extraction. Sampling was carried out under permits from the Department of Agriculture, Forestry and Fisheries (RES2010/01, RES2011/12).

DNA extraction, PCR amplification and sequencing
Genomic DNA was extracted from cirral tissue using the Wizard SV genomic DNA purification system (Promega). PCR amplification and sequencing of a fragment of the COI gene was initially carried out using universal primers [49]. The data generated by these primers were used to design species-specific degenerate primers which were subsequently used in all reactions (COI-TsF: 59-CGG RGC TTG ATC CGC YAT RGT MGG-39 and COI-TsR: 59-GCT CCG GCT AAA ACT GGA AKH G-39). Nuclear DNA data were generated from a subset of samples by amplifying the ITS1 region (including partial sequences of the 18S and 5.8S genes that flank the ITS1 region) using the primers of Bass et al. [50]. Amplification was performed in 50-ml reactions containing 4 ml of 10-100 ng DNA, 2.5 ml each primer (0.1 mM), 5 ml 10xPCR reaction buffer (Super-Therm), 5 ml dNTPs (0.2 mM, Thermo Scientific), 4 ml Mg 2+ (Super-Therm), 0.2 ml Taq polymerase (Super-Therm) and 13.4 ml distilled water. PCR reactions were carried out in an ABI GeneAmp 2700 automated thermocycler with the following conditions: an initial denaturation step at 94uC for 3 min, followed by 35 cycles of denaturing at 94uC for 30 s, annealing (for 30 s) and an extension at 72uC for 45 s, ending with a final extension at 72uC for 5 min. PCR annealing temperatures were set at 58uC and 56uC for the COI and ITS1 loci, respectively. The amplified products were gel separated and purified using the Wizard SV Gel PCR Clean-Up System (Promega). Purified products were cycle-sequenced using BigDye (Applied Biosystems [ABI]) and analyzed on a 3730 automated sequencer (ABI).

Genetic diversity
Sequences were manually inspected and edited using BioEdit v7.0.9.1 [51], and the mtDNA nucleotides were translated into protein codons to confirm functionality. All haplotypes have been deposited in GenBank with the following accession numbers KC935448-KC935857 (COI) and KC967663-KC967682 (ITS1).
Haplotype designation of the nuclear DNA data was performed in PHASE, using a coalescent-based Bayesian method [52] as implemented in DnaSP v5.1 [53]. A total of 1 000 000 iterations were performed using default settings. Arlequin 3.5 [54] was used to calculate the number of polymorphic sites and the haplotype (h) and nucleotide (p) diversities for both marker types. [54]. Differentiation across sampling localities was examined through pairwise W ST comparisons. Bayesian geographical structuring was investigated using BAPS v5.4 [55], where the geographic coordinates of each sampling locality were obtained from Google Earth. To investigate geographic distance as a contributor towards genetic structure, isolation by distance [56] was tested using a Mantel test [57] as implemented in Arlequin 3.5 with 1000 permutations. The geographical distances between sampling localities were measured as the shortest land-sea interface in Google Earth.
Statistical parsimony haplotype networks were constructed using TCS 1.21 [58] and haplotypes were connected using a 95% confidence interval. The higher level clustering of clades that were not connected in the haplotype networks was depicted using Bayesian analyses, coupled with the Markov Chain Monte Carlo (BMCMC) inference, conducted in MrBayes v3.04b [59]. The most likely (best probabilistic) model of sequence was determined using the Bayesian Information Criteria (BIC) and Akaike Information Criteria (AIC) for the COI and ITS1 gene, respectively, as employed in the program jModelTest v0.1.1 [60]. Bayesian analyses of the COI and ITS1 sequence data were run for 10 million and five million generations, respectively, with a burn-in of 25%. For each gene, two different runs were conducted to test for convergence, each with four chains and sampling frequency of 1000. Thereafter, the post-burn in samples of the two runs were combined and a 50% majority-rule consensus tree was constructed in MrBayes v3.04b. Clade support was assessed using posterior probabilities for the Bayesian tree and only values $0.95 were regarded as robust support. Tetraclita squamosa and Tetraclita pacifica were used as outgroup taxa (GenBank accession numbers DQ647770 and DQ363746, and DQ363680 and DQ363740, respectively).

Demographic history
Arlequin 3.5 was used to test the influence of historical sea level changes on the demography of T. serrata. Fu's F s test [61] was performed to test for deviations from neutrality. To confirm that departures from neutrality were caused by demography and not by natural selection, the McDonald-Kreitman (MK) test was performed in DnaSP v5.1, using T. squamosa and T. pacifica (GenBank accession numbers DQ363697-DQ363706 and DQ363686-DQ363695) as outgroups, respectively. This was followed by a mismatch distribution to test for population expansion and to estimate parameters needed to date possible  Table 1) for Tetraclita serrata are shown in addition to coastal regions [west (W), southwest (SW), south (S), east (E)] in which samples were collected. These coastal regions correspond broadly with the Namaqua, south western Cape, Agulhas and subtropical Natal Bioregions, respectively [47]. The frequency of Clade A and Clade B within each sampling locality are indicated using pie charts. The major ocean currents are shown. The direction and intensity of gene flow between sampled localities, and the overall direction of gene flow are indicated by arrows; mutation-scaled migration rates between localities are indicated above these arrows. doi:10.1371/journal.pone.0102115.g001 population expansions [54], [62]. The formula T = t/2 mk [62] was used for dating approximations. As the exact mutation rates and generation times for T. serrata are unknown, we used the generation times of various barnacle species, ranging from one to two years [8], [32], [63], [64], and two published COI lineage mutation rates for barnacles (1.55610 28 [65] and 2.76610 28 substitutions per site per generation [18]). A Bayesian Skyline Plot (BSP) [66] was constructed using BEAST v1.8 [67] to estimate changes in population size over time and to date a possible population expansion. The unique haplotypes of each COI data set (Clades A and B) were run for 300 million and 200 million generations, respectively, under a GTR+I+G nucleotide substitution model with individual parameters estimated from the data and a constant skyline model with five groups. For each BSP, prior distributions for the root height of the population were notified by initial estimates from the mutation rates used in the mismatch distribution analyses. The chain was sampled every 100 000 generations and the first 10% was discarded as burn-in. Trace plots were inspected to assess convergence, mixing and stationarity of the MCMC process in Tracer v1.6 [68]. The effective sample sizes were checked and confirmed to be $200 in order to avoid autocorrelation of parameter sampling.

Connectivity and direction of relative gene flow between sampling localities
The program Migrate 3.1.6 [69] was used to calculate the relative migration rates between sampling localities inferred from the mtDNA data. Given the linear South African coastline, a stepping-stone model was used [70]. The settings were as follows: 10 short-chains, each with a total of 2,500 generations and a sampling increment of 20 generations, and 3 long-chains, each with a total of 50,000 generations and a sampling increment of 50 generations. The first 10,000 genealogies were discarded (burn-in). An adaptive heating scheme with four chains (starting values of 1.00, 1.50, 3.00 and 6.00) and a swapping interval of one was used to ensure that efficient mixing occurred. For the other settings, default values were implemented. The analyses were run twice to ensure congruence.

Genetic diversity and population structure
Mitochondrial DNA analyses. The COI data analyses were based on a final alignment of 499 bp for 410 individuals. The maximum parsimony network, Bayesian tree and spatial clustering of all individuals (BAPS) indicate that T. serrata comprises two genetically distinct groups with strong nodal support. The two genetic assemblages occur sympatrically (Fig. 2a-c). These haplogroups represent two reciprocally monophyletic lineages, both showing Bayesian Posterior Probabilities of 1.00 (BPP). These are labeled Clades A and B (Fig. 2) and are separated by a mtDNA sequence divergence of ,8%. In proportion to Clade A, the frequency of individuals belonging to Clade B is markedly smaller (a sampling ratio of 12:1, Fig.1).
Within Clade A, haplotype diversity was 0.979 (60.003), with 187 haplotypes detected, and within Clade B it was 0.976 (60.019), with 26 haplotypes detected. The nucleotide diversity of the complete data set was 0.019 (60.01); generally mtDNA genetic diversity was highest on the south-west and south coasts, with a progressive decrease towards the east coast and the lowest values on the west coast (Fig. 3a). With the removal of Clade B, the nucleotide diversity values are lower, with the south-west coastal region generally being the highest (Fig. 3). Although the haplotype diversities determined for each locality for the entire data set (Clade A and Clade B individuals) were generally high (.0.94, Fig. 3b), on a broader scale, sampled localities on the south-west coast, as well as the closest localities sampled on either side of the south-west coast region, namely Jacob's Bay, Kommetjie and Cape Agulhas, each had more than 50% unique haplotypes (Table 1). In contrast, those haplotypes shared among regions occur more frequently along the east coast than in the other regions (i.e. fewer private haplotypes, less than 45%).
Within Clade A, the haplotype network (Fig. 2b) shows that a considerable amount of haplotype sharing is found along the east coast. Also, west and south-west coast sampled localities group together and are more divergent, suggesting spatial genetic structure. In Clade B, the majority of haplotypes are representative of individuals sampled along the west, south-west and south coasts, with only two Clade B individuals sampled as far as Haga Haga. Many of the haplotypes within this haplogroup are separated from the common haplotype by an intermediate ''unsampled or extinct'' haplotype (Fig. 2b).
AMOVA suggests shallow, but significant, genetic differentiation among the 20 sampled localities (W ST = 0.035, P,0.01); most of the variation was present within sampled localities (96.44% for the entire data set). When sampling localities were partitioned into their respective coastal regions (Table 1), AMOVA revealed significant differentiation among groups (W CT = 0.041, P,0.01). When the individuals belonging to Clade B (Fig. 2b) were removed from subsequent analyses, the W ST value increased (0.065, P, 0.01). For the entire data set, pairwise W ST analysis showed that Kommetjie was significantly differentiated from all localities east of Cape Agulhas (south coast), with W ST values ranging from 0.09 to 0.23 (south coast, Table 2). Kommetjie was also found to be significantly differentiated from most of the localities sampled on the west coast (Table 2). Betty's Bay was found to be significantly differentiated from all localities east of Herold's Bay (south coast), with W ST values ranging from 0.1 to 0.19 (Table 2). Within Clade A, Kommetjie and Betty's Bay were significantly differentiated from all sampling localities, except from each other and from Port Nolloth, with W ST values ranging from 0.07 to 0.34 (Table 2). Due to the relatively small sample size represented by individuals of Clade B, pairwise W ST analysis was not conducted for this clade.
Nuclear DNA analyses. The aligned sequence data for ITS1 were 426 bp in length, with 31 unique alleles and four indels detected. The ITS1 data displayed a high haplotypic diversity of 0.960 (60.019) and a low nucleotide diversity of 0.010 (60.005), which is comparable with that of the mtDNA data. No significant structure was obtained and Clades A and B formed by the mitochondrial DNA analyses collapsed into a single haplogroup (Fig 2d, e).

Isolation by distance
A significant correlation between geographic distances and genetic differentiation was obtained for the entire T. serrata data set (Mantel test, Z = 8.925, P,0.01), revealing a larger Z-value of 14.176 (P.0.05), which was not significant, with the removal of Clade B. The Mantel test was found to be non-significant when isolation by distance (IBD) was tested within each coastal region [west coast, Z = 0.119 (P.0.05); south coast, Z = 20.035 (P. 0.05); east coast, Z = 20.096 (P.0.05)]. South-west coast localities were not included in the test because only two localities had been sampled in that region.

Migration estimates between sampling localities
Within the sampling range of T. serrata, the coalescent stepping-stone model revealed migration rates that were largely bi-directional in the south and east coastal regions, whereas a more asymmetrical gene flow pattern was evident in a westerly direction between Mossel Bay and Hondeklip Bay (with the exception between Betty's Bay and Cape Agulhas) (Fig. 1). Interestingly, no southward gene flow was detected between Jacob's Bay and Betty's Bay, around Cape Point (Fig. 1).

Demographic history
For the mtDNA data, Fu's F s revealed negative values for the entire data set (although it was not significant for all sampling localities; Table 1), Clade A (F s = 217.858, P,0.01) and Clade B (F s = 225.065, P,0.01), respectively. The MK test was not significant for T. serrata when T. squamosa (G-value = 0.841, P.0.05) was used as an outgroup. The more closely-related species T. pacifica could not be used for this test because of the lack of nonsynonymous fixed differences with T. serrata in the  Table 1), which tests for a fit to a unimodal mismatch distribution.
The time since population expansion was calculated for mtDNA Clades A and B ( Table 3). Analyses of Clade A suggest the time since expansion between ,40 000 and 145 000 years ago. Similarly, for Clade B, the time since expansion was calculated to be between ,46 000 and 164 000 years ago. Bayesian Skyline Plot analyses reveal dates similar to the mismatch distribution analyses. Analyses of Clade A show a population expansion took place between ,40 000 and 100 000 years ago, while analyses of Clade B show a population expansion took place between ,50 000 and 100 000 years ago (Fig. 4).

Within-population diversity
Genetic diversity indices varied across the sampling range, with striking differences between sampling areas; generally, diversities were highest on the west and south-west coasts, and they decreased towards the eastern localities. A distinct peak in diversity was found for the south-west coast (Fig. 3). There were also notable differences in the distribution of private/unique haplotypes, with a higher number of private haplotypes found in the south-west coastal region ( Table 1).
Regions of higher genetic diversity can be explained by several scenarios. Genetic variation can accumulate over evolutionary time if a region contains an area of high local retention or larval accumulation, especially where larvae accumulate in the nearshore environment [29], [71], [72]. High genetic diversity is also expected in a region where a large population has been maintained over evolutionary time, possibly due to the stability of the habitat [73]. If contemporary population size is an indicator of evolutionary history, then this will have contributed to the high genetic diversity as T. serrata usually maintains large population sizes on most rocky shores [45], [74], [75]. Irrespective of population size, regions of higher genetic diversity can also signal refugial areas during glacial low sea-level stands [76]. It can also be a signal of the absence of strong directional or stabilizing selection, as both mechanisms have been shown to reduce genetic diversity [77]. Populations at or near refugia should retain genetic diversity and may also be indicated by the presence of ancestral haplotypes [78]. Interestingly, Betty's Bay on south-west coast, possesses the highest combination of haplotype and nucleotide diversity (Tables 1, 2). A study on the commercially exploited rock lobster, Jasus lalandii, also identified the region between Betty's Bay and Kommetjie as one of high diversity [21], but the processes driving this area to be more genetically diverse than others are not understood.

Genetic structure of Tetraclita serrata
Population genetic structure of T. serrata is characterized by a significant AMOVA and a strong signal of isolation by distance for the mtDNA data. An important result was that T. serrata is comprised of two genetically distinct groups that are geographically mixed (Fig. 2b, c), consistent with Tsang et al. [46]. Although Clade B is highly divergent and shares no haplotypes with Clade A at the mitochondrial level, the shallow population structure recovered reflects the geographic overlap of the haplotypes of the two clades, with the majority of the haplotypes being shared among the sampled areas (Fig. 2b). The pattern observed reflects a significant historical differentiation of the two clades and may represent an incipient speciation event. Lowered sea levels, as a consequence of historic glaciations, may have been responsible for the formation of allopatric lineages that have since come into secondary contact.
In the absence of a clear allopatric, geographic explanation that could have given rise to the two clades in T. serrata, it is important to consider that physiological differences could also result in genetic diversification. Field observations and physiological experiments on the mussel Perna perna indicated that the western lineage of the species is intolerant of high water temperatures, thus explaining its exclusion from the east coast [79]. Similarly, in the mudprawn Upogebia africana, larvae of the eastern lineage are unable to tolerate cooler waters and, thus, are unable to establish themselves in cool-water provinces [80]. The mtDNA data of T. serrata suggest that Clade B extends no further east than the border of the sub-tropical and temperate regions at Haga Haga, whereas Clade A can be found in the sub-tropical region (Fig. 1,  Fig. 2b). This implies that Clade B may be constrained to more temperate environments, as a result of physiological intolerance. The annual mean sea surface temperatures experienced along the South African coastline ranges between 12uC on the west coast to 26uC on the east coast [81]. Tetraclita serrata displays regular and annual recruitment and Dye [82] revealed a high variation in spatial and inter-annual recruitment intensity along the east coast of South Africa. It was also found that recruitment density was related to mean abundance [82]. This suggests that larval dispersal  Table 2. Pairwise WST values among sampling localities for T. serrata using the COI gene.  and settlement may differ among regions. Along these lines, it is particularly notable that two distinct, species-level lineages of the Red Sea barnacle Tetraclita rufotincta occupy different vertical zones of the intertidal [83]. The genetic divergence in T. serrata may be explained by selective pressures due to ecological conditions, such as desiccation tolerance and predation. Given the geographic overlap of the two clades, it is possible that the genetic pattern found in T. serrata is the result of an ecological speciation process. Individuals belonging to Clades A and B may occupy different zones of the mid-intertidal (micro-habitat isolation), although this remains to be tested further.
Despite the clade structure present in the mtDNA data and a sequence divergence of ,7.92%, the lack of geographic structure for Clades A and B and the overall lack of structure for the nuclear data suggest incipient speciation at best (Fig. 1, 2). This result is supported by Tsang et al. [46], who similarly found no genetic  structure for nuclear DNA data (Histone 3) and no morphological differentiation. As the species is hermaphroditic, it can be ruled out that the lack of nuclear genetic structure is that of sex-biased dispersal, and is more than likely due to differences in the fixation rates of the two marker systems.

Population genetic structure: associations with lifehistory characteristics
Several phylogeographic studies for southern African marine species suggest the presence of at least four marine phylogeographic breaks (see [19], [25]), of which we sampled across three (Cape Point, Cape Agulhas and the area between Port Elizabeth and the northern Transkei; see Fig. 1 of [19]). Although globally many studies have shown that life history characteristics such as pelagic larval duration (PLD) are not correlated with population genetic structure and gene flow [84], [85], planktotrophic larval dispersal of T. serrata may give some insight into the lack of structure recovered in this study. Even though the shared haplotypes found for the entire data set indicate that this brooding barnacle is able to disperse over long distances, suggestive of a long PLD, T. serrata does not form a panmictic population as several sampled localities are differentiated, even if W ST values are low ( Table 2). The presence of large numbers of haplotypes unique to certain regions (Table 1), particularly in the south-west coastal area (KO, WP, BB), provides evidence for population genetic structuring in the west and south-west coastal regions. The significant differentiation of Kommetjie (Cape Peninsula) and Betty's Bay (False Bay region) from other localities (Fig. 2b, Table 2) suggests that, at the finer geographic scale, dispersal may be restricted by water currents or by the local retention of larvae in the water column, as has been shown for other species [34], [86], [87], [88]. This is further supported by the signal of IBD for T. serrata.

Gene flow in relation to ocean currents
The strength and directionality of gene flow in T. serrata show interesting correlations with oceanographic features along the South African coast. All west coast (except Port Nolloth) sampled localities experience gene flow primarily in the direction of the Benguela Current, while half of the localities sampled on the south-west coast (see below) experience gene flow in the direction of the Agulhas Current (Fig. 1). Ocean currents and prevailing wind direction have long been recognized as important factors affecting the distribution, abundance and genetic variation of marine benthic invertebrates with planktonic larvae [89], [90]. A pattern of unidirectional gene flow, which is much lower than relative gene flow estimates around Cape Point, may play a role in the genetic structuring of sampled localities on the south-west coast, evident in the haplotype network with haplotypes from the west coast and south-west coast sampling localities clustering together (Fig. 1). Such reduced gene flow was also shown for rocky shore fishes [22], [26] and an endemic sea urchin, Parechinus angulosus [24] strongly suggesting that gene flow is reduced for inshore coastal species in this region. Cape Point is also characterized as a region where a major biogeographic break occurs [19]. The barrier is thought to be the result of the distinct temperature regimes on either side of the peninsula, which may influence larval dispersal and/or survival in the different environments. By contrast, the signal is more mixed on the south and east coast, with bi-directional gene flow, suggesting that both the inshore, wind-driven countercurrents and the offshore Agulhas Current play a role in the distribution of larval T. serrata [19].

Demographic history
The combination of high haplotype diversity and low nucleotide diversity is a typical signature of population expansion following a bottleneck [91], [92]. Fu's F s simulations (Table 1) and the MK test, which was not significant, indicate that the observed deviations from equilibrium are not caused by natural selection and that T. serrata has undergone recent demographic change. The mtDNA haplotype network, showing several common haplotypes (Fig. 2), mismatch distributions (Table 1) and the timing of the most recent expansion event for T. serrata are indicative of a growing population. The population expansion of T. serrata predates the LGM, taking place 164 000 to 31 000 years ago, depending on the mutation rate and generation time used ( Table 3). The BSP analyses agree with this (40 000 to 100 000 years ago, Fig. 4). Although many southern African marine taxa show patterns of post-LGM population growth [8], [20], [21], [93], [94], [95], [96], also confirmed by a number of barnacle species found elsewhere [5], [8], several studies among rocky shore marine species also show pre-LGM population expansions (Pacific Ocean: [8], [18], [97], [98], [99], [100]; Atlantic Ocean: [4]). This result suggests that the dramatic and frequent climatic changes of the Pleistocene, pre-dating the LGM, may have had a stronger impact in shaping the demographic and biogeographic history of T. serrata [8], [101]. Even if a more conservative approach was used, such as the mitochondrial mutation rate based on the divergence rate of 2% per million years [102], the demographic changes of T. serrata would still have predated the LGM. The two mutation rates that were used from barnacle genera were substantially greater than the 2% rate and still the population expansion of T. serrata predated the LGM. Therefore, despite the limitations of calculating the timing of expansion [27], the variation in mutation rates alone cannot account for the relatively deep demographic history for T. serrata.
Broadly, our results support Marko et al. [8] in that suspension feeders possessing planktonic larvae that occur in the midintertidal region show patterns of LGM persistence. The majority of studies to date on South African marine species showing patterns of post-LGM recolonization also possesses planktonic larval stages and is found subtidally, and is neither suspension feeders nor sessile species [7]. Nevertheless, these are certain traits worth considering for future research given that there are a number of barnacle species that possess planktonic larvae and live in rocky shore habitats showing clear population genetic signatures of regional persistence [4], [8], [18], [98].