Fine ‐ scale biogeography: tidal elevation strongly affects popula ‐ tion genetic structure and demographic history in intertidal fishes

. Numerous studies have demonstrated population genetic structuring in marine species, yet few have investigated the effect of vertical zonation on gene flow and population structure. Here we use three sympatric, closely related clinid species, Clinus cottoides , C. superciliosus and Muraenoclinus dor ‐ salis , to test whether zonation on South African intertidal rocky shores affects phylogeographic patterns. We show that the high ‐ shore restricted species has reduced gene flow and considerably higher F ST val ‐ ues ( F ST = 0.9) than the mid ‐ and low ‐ shore species ( F ST < 0.14). Additionally, we provide evidence for remarkably different demographic and evolutionary histories, ranging from extreme population bottle ‐ necks to population persistence, which are probably linked to effective population size and habitat spe ‐ cialisation. This study further highlights the need for a multispecies approach to unravel the biological and evolutionary processes that drive extant population genetic patterns in marine species, as even closely related species with similar life histories show highly variable results.


Introduction
Marine rocky shores encompass some of the most dynamic environments in the world, with animals and plants having to survive extremely variable and demanding conditions (Denny and Gaines 2007). Fauna and flora of all rocky shores show a degree of vertical zonation, extending from the subtidal to the upper margins of the shoreline. Zonation depends on several complex physical and biological interactions, including elevation of the shore and wave exposure as well as the physiological capacity of intertidal organisms (Dayton 1971, Blamey andBranch 2009).
Molecular tools have greatly enhanced our understanding of the population structuring of many rocky shore inhabitants, as well as helped to elucidate the evolutionary history of many marine species (e.g., Ayre et al. 2009, Pelc et al. 2009, von der Heyden 2009). However, one question that remains poorly understood is whether sympatrically occurring rocky shore species, which are predominantly found at different shore heights, exhibit differences in population genetic structure across their distributional range. Many studies have focussed on the duration of pelagic larval dispersal in relation to population genetic structuring of marine species (e.g., Waples 1987, Doherty et al. 1995, Selkoe and Toonen 2010, Riginos et al. 2011), but whether vertical zonation of marine species on rocky shores also results in different population genetic structures remains uncertain. tidal elevation affects population genetic structure Theoretically, species restricted to the high shore may show greater levels of population structuring, as the potential for dispersal may be more limited compared to organisms that inhabit the mid to low shore (all other things being equal). This is a function of time spent submerged during tidal cycles, and thus time available for movement, spawning and the subsequent broadcasting of propagules. In a study focussing on 50 intertidal species, Kelly and Palumbi (2010) suggest that species in the mid to high intertidal zone are more strongly genetically structured than other rocky shore taxa. In New Zealand triplefins, population structuring also appears to be correlated with depth of habitat (Hickey et al. 2009) and in two species of mudsuckers (Gillichthys spp.), the highshore G. seta showed less gene flow than did G. mirabilis, which occurs lower in the intertidal zone (Huang and Bernardi 2001). In contrast, Bird et al. (2007) demonstrated that in three Hawaiian Cellana limpets, C. talcosa, which is restricted to the shallow subtidal, showed the strongest population structuring.
Further, molecular tools have greatly facilitated the reconstruction of estimates of quantitative demographic change of species distribution and abundance, although this has not been without pitfalls (see Karl et al. 2012). There is evidence from southern African marine species that paleoclimatic changes have played significant roles in shaping their demographic and population genetic patterns, with many species showing signals of population contraction and expansion within the past 100,000 years (von der Heyden et al. 2010, Teske et al. 2011, Muller et al. 2012. While two sympatric limpets of the genus Siphonaria show similar phylogeographic breaks (Teske et al. 2012), the processes that led to the formation of the breaks were quite different. Therefore, unravelling past demographic histories is important for explaining extant population patterns, and phylogenetically closely related species with sympatric distributions are ideal models for testing the effects of past climatic change.
The family Clinidae has a globally disjunct distribution (von der Heyden et al. 2008). The cool -temperate environments of the western and southern South African coasts are the hot-spot of clinid diversity, with at least 50% of all clinid species and several cryptic species and species complexes occurring there , Holleman et al. 2012. As is the case for most specialised intertidal fishes, clinid fishes are bottom-dwelling and cannot sustain prolonged periods of swimming. Dispersal ability in South African clinid fishes is probably limited, as all are viviparous, giving 'birth' to well-developed postflexion larvae. Vivipary should theoretically limit the dispersal potential of clinid fishes, especially as it is unlikely that adult fishes move great distances. As expected, the South African clinids, Clinus cottoides (Valenciennes, 1836), C. superciliosus (Linnaeus, 1758) and M. dorsalis (Bleeker, 1860), show strong population structure along ~1,500 km of the southern African coastline (von der Heyden et al. 2008. This study further explores phylogeographic structuring of three clinid species, C. superciliosus, C. cottoides and M. dorsalis, which differ in their vertical distribution as mature adults on South African rocky shores. Clinus superciliosus [total length (TL) = 30 cm] is the most abundant clinid on rocky shores, where adults are much more abundant on the lower shore, with juveniles occurring primarily on the higher shore (Prochazka and Griffiths 1992). Clinus cottoides (TL = 15 cm) is more abundant in the mid-shore area than either the low or the high shore (Bennett and Griffiths 1984). Muraenoclinus dorsalis (TL = 10 cm) is a remarkable fish in that it can spend time out of water at low tide, with individuals often sheltering under damp seaweed or rocks (authors pers. obs.). Mature M. dorsalis are found predominantly in the 'upper balanoid zone', which lies farther up the vertical shoreline than the areas where adults of the other two species are found. Therefore, reproductively active adults of all three species show distinct vertical zonation on South African shores (Prochazka 1996).
This study has three major aims: 1) to determine whether the high shore dwelling M. dorsalis exhibit stronger population genetic structuring than the mid or lower shore species, 2) to determine whether patterns of gene flow among intras-pecific populations for each of the three species of clinids are congruent and 3) to determine whether each of the species shows similar patterns of demographic history. Based on theoretical expectations, we hypothesise that greater genetic structuring should be apparent for the high-shore M. dorsalis than for the mid-or lower-shore inhabiting C. superciliosus or C. cottoides, respectively, with gene flow in M. dorsalis being much more restricted than for the other two species. Finally, we hypothesise that each of the species will show signals of demographic expansion following a population bottleneck, as this pattern is com-monly recovered for a variety of other southern African marine fish species and may well be linked to climate driven changes during recent glacial cycles (von der Heyden et al. 2010).

Sampling and molecular protocols
Fishes were sampled throughout the core of their distributional ranges in the south-eastern Atlantic between Jacobsbaai on the west coast and Gansbaai on the south-west coast of South Africa   Table 1. Coalescent estimates of migration rates for Nm are also given between adjacent populations as calculated from a stepping stone model. Numbers above the arrow give Nm from east to west or south to north, whereas the number below the arrow is from west to east or north to south. respectively, were sampled (Table 1); sequences for C. cottoides were taken from von der Heyden et al. (2008). Fishes were collected exclusively from tidepools using hand nets and the anaesthetic Ethyleneglycolmonophenylether (Merck). Muscle tissue was dissected from each fish and total genomic DNA extracted with the DNEasy kit (Qiagen) following the manufacturer's protocol.
PCR primers and conditions to amplify the variable 5' region of the mtDNA control region were as described by von der Heyden et al. (2008). PCR products were gel purified using the GfX kit (Amersham Biosciences). Purified products were sequenced using Big Dye Terminator chemistry (Applied Biosystems) and run on an ABI 3100 automated sequencer. All sequences were submitted to GenBank with the following accession numbers KC577471-KCKC577492 (M. dorsalis), KC577493-KC577538 (C. superciliosus).

Population structuring and gene flow
Sequences were edited and aligned in BioEdit (Hall 1999). For each species, data were analysed in two ways; all individuals in one data set (i.e. all sampling localities pooled for each species), as well as three data sets in which individuals for each species were divided into sampling localities. Standard diversity indices [haplotype (h) and nucleotide (π) diversity] were calculated for each data set. The lowest and highest effective number of haplotypes was also calculated for each species using the approach of Jost (2007). To investigate population genetic structuring amongst sampling localities, analysis of molecular variance (AMOVA) and pairwise population Φ ST values were carried out in Arlequin 3.5 (Excoffier and Lischer 2010). Significance levels were calculated with a nonparametric permutation approach with 10,000 iterations. To visualise the relationships among haplotypes, we made use of a statistical parsimony network as implemented in TCS 1.21 (Clement et al. 2000).
Gene flow is an important contributor to population genetic patterns. Therefore we calculated a mean value of gene flow, Nm, between sampling localities using a stepping stone model  Felsenstein 1999, 2001). We used the settings as in von der Heyden et al. (2008), but ran each analysis three times to check for consistency and averaged the results across runs. We also used Gelman's convergence criterion to check for convergence between chains.

Historical demography
To investigate the demographic history of the three species, Fu's F S was used to test for demographic change. These tests were carried out in Arlequin 3.5. To decouple the potential for selection and historical demography, we evaluated population fluctuations based on coalescent models for each species separated into five or six populations. Population characteristics and sample numbers are given in Table 1. Population parameters Θ (theta) = 2Nμ, where μ is the mutation rate for mtDNA and g (the exponential growth parameter in units of μ) were estimated using a coalescent approach with Fluctuate 1.4 (Kuhner et al. 1998). The parameters Θ and g were estimated jointly. Estimates were obtained by running four replicates, which generated a mean value and its associated standard deviation. Analysis of each data set was conducted with 10 short Monte Carlo chains of 4,000 steps each and five long chains of length 20,000, with a sampling increment of 20. Fluctuate generated a random topology for initial searching. Population coalescence times were also estimated by assuming that coalescence was reached when the population size was reduced to 1% of its present-day value, following Wares and Cunningham (2001). In order to estimate coalescence time in the absence of a molecular clock calibration for South African clinids, we used the mutation rate (μ) as μ = substitutions per site per generation obtained for trans-isthmian geminate fish species in the genus, Anisotremus (Bernardi and Lape 2005), which range from 1.7x10 -8 to 0.9x10 -8 mutations per site per generation. An AMOVA revealed different levels of genetic structuring for the three species; M. dorsalis had the highest global Φ ST (Φ ST = 0.9, P < 0.0001), whereas the Φ ST for the two Clinus species was similarly low (Φ ST C. cottoides = 0.05, P < 0.0001; Φ ST C. superciliosus = 0.07, P < 0.0001). Significant pairwise Φ ST values ranged between 0.17-0.96 for M. dorsalis, between 0.07-0.14 for C. cottoides and between 0.08-0.11 for C. superciliosus ( Table  2). The haplotype networks clearly showed the relationships amongst fish sampled at different localities ( Figure 2). Haplotype networks for the two Clinus species predominantly show shared haplotypes among different sampling localities. In contrast, the network for M. dorsalis showed mostly haplotypes that were unique to specific locales, a result that is consistent with the high genetic differentiation recovered. Evaluation of the limits of statistical parsimony suggests that topologies connecting haplotypes by eight steps or fewer have a cumulative probability of 95% of being correct. For M. dorsalis, however, not all haplotypes could be connected under a 95% confidence limit.

Estimation of gene flow
Coalescent analyses suggested that gene flow, Nm, was limited among populations for all three species (Figure 1). Again, the most restricted levels of gene flow were recovered between M. dorsalis populations, whereas over the same spatial distances, values between C. cottoides popula-  1-3) was strongly asymmetrical, with northward but not southward dispersal. On the south-west coast (localities 4-6), gene flow was more bi-directional, although generally the values of Nm were low (Figure 1).

Estimates of population growth
Estimates of both population size (Θ) and population growth (g) (

Comparison of genetic structuring and gene flow in live-bearing intertidal fishes
This study compared three species of clinid fishes across their core distributions in south-western Africa. As the females of each species give birth to live young that probably recruit into natal rock pools, we can discount larval dispersal as a factor shaping the genetic structure of each fish species. This provides us with an opportunity to investigate the relative importances of other factors that may explain phylogeographic and demographic patterns in this dynamic marine environment.
Our first hypothesis stated that we expected the high-shore dwelling M. dorsalis to display higher levels of population genetic structuring compared to the two Clinus species that are abundant from the mid intertidal to the shallow subtidal. This is indeed supported by our analyses, with M. dorsalis displaying high pairwise Φ ST values and a global F ST of 0.90, compared to the lower values of C. cottoides and C. superciliosus (Table 2a,  Gansbaai -86.9 ± 2.9 0.005 ± 0.00 2305 ± 1779 0.008 ± 0.006 178.0 ± 35.6 0.05 ± 0.024 Sampling locality Table 3. Historical demography of South African clinids. Left column gives sampling locales. Each locale has two rows: the first is growth (g) the second row is theta (Θ) calculated with variable growth. Θ values correspond to averages of four replicates and its standard deviation.
immediately adjacent populations (Figure 2). It is likely that the high population genetic structuring is driven by the ecology of this high-shore species; M. dorsalis is able to use damp rocks and seaweeds on the high-shore, where no other rocky shore fishes are able to compete (Prochazka and Griffiths 1992). Per our expectations, this lifehistory strategy greatly reduces the dispersal ability of M. dorsalis, which experiences only tidal influx for a few hours each day, compared to adult C. superciliosus, which are found in the shallow subtidal environment and therefore have a much greater potential to disperse. A similar pattern was described by Kelly and Palumbi (2010): in their meta-analysis high-and mid-intertidal species had larger F ST values compared with species on the lower shore or in deeper water. Body size and environmental tolerance are also important in allowing fishes to move across barriers; larger fishes with broad environmental tolerances are more likely to cross oceanic barriers than smaller or more specialised fishes (Bradbury et al. 2008, Luiz et al. 2012. Although little is known about the ecology of clinid fishes, ecology may be an important consideration, as M. dorsalis is the smallest of the three species. No differences were apparent in the genetic structures of the mid-tidal and low-tidal Clinus species, with both species showing similar measures of genetic differentiation (Table 2). This is somewhat surprising as adult C. cottoides are seldom found in sub-tidal or low shore habitats and should therefore be more movement restricted than C. superciliosus. However, Clinus cottoides has the lowest genetic diversity of the three species, with over 72% of individuals sharing a common haplotype. This genetic signature may be indicative of extreme population reduction in the past, with subsequent recolonization of parts of its range, which limits the inferences that can be made from a single mtDNA locus.
Gene flow was generally low among populations of all three species, with broadly congruent patterns. For example, there is limited-to-no southward gene flow from Jacobsbaai to Kommetjie along the west coast, but some in the opposite direction. Along the south-west coast (Wooley's Pool to Gansbaai), gene flow is bi-directional, with no distinct pattern. Notably, M. dorsalis has the smallest amount of gene flow between populations of the three species, leading to high levels of population differentiation. Mid-shore Clinus cottoides had the highest levels of gene flow, which is surprising as we expected C. superciliosus to show higher levels. Again, this is probably affected by the strong population bottleneck C. cottoides experienced (for more discussion see below).
From a broader biogeographic perspective, it appears that none of the species shows particularly strong population structure across one of the most pronounced biogeographic breaks, Cape Point, on the south-west coast of Africa (Griffiths et al. 2010). There is some evidence for other South African marine species that biogeographic and phylogeographic breaks are congruent (see Teske et al. 2011), but this is not the case here, as none of the species shows its largest F ST between Kommetjie and Wooley's Pool or Betty's Bay. Rather, for M. dorsalis, a more pronounced break appears to lie between Kommetjie and Sea Point (a distance of ~27 km), whereas for C. cottoides, the greatest genetic differentiation lies farther north between Jacobsbaai and Sea Point.

Differing evolutionary patterns in three sympatric rocky shore fishes
Our third hypothesis, that co-distributed and closely related species show similar patterns in their demographic history, is not supported. First, our results show that the three clinid species have different signals of demographic history. Midshore C. cottoides has much lower haplotype and nucleotide diversities than the other specides (Table 1), and despite the greater number of fish sampled, it had only 16 haplotypes. This, in conjunction with the star-burst haplotype network and the largest signal of growth, suggests that C. cottoides likely underwent an extreme bottleneck in population size, followed by demographic expansion. In contrast, the haplotype network for C. superciliosus (Figure 2) suggests that this species has a much longer evolutionary history, given the number of older haplotypes that are characterised by several mutational step differences, which are absent in the C. cottoides network (Figure 2). This is also particularly evident in coalescent estimates of growth and theta, which showed that C. superciliosus had far longer coalescence times than the other two species.
Second, this pattern might be explained by environmental tolerance; C. superciliosus, which is a generalist and thus able to use a greater portion of available habitat, may have adapted better to past drivers of change, whereas the mid-shore restricted C. cottoides may have been more susceptible to environmental change. Clinus superciliosus is also twice as abundant in intertidal pools as C. cottoides or M. dorsalis (Bennett andGriffiths 1984, Prochazka 1996), which means that new mutations will be retained for longer periods of time. This assumes that the effective population size of C. superciliosus is also larger than that in the other two species, and values of theta (Table 3) support this scenario.
Third, shallow subtidal habitat is more continuous than rocky shore, so species that rely (i.e. are more specialised) on rocky shores may have had populations extirpated when sea-level dropped during glacial cycles. There might have been significant rocky shore habitat loss, especially on the south coast (J. Compton, UCT, pers. comm.), and this could have led to population decreases in rocky shore specialists. This may explain the pattern for M. dorsalis, which is the most specialised of the three species and which does not show a signal of demographic change in neutrality tests (Table 1). Coalescent results are more mixed, but indicate that at least some populations experienced low or negative growth. Given the extreme high-shore environment M. dorsalis at present occupies, this species may well have a more specialized physiology that facilitates adaptation to anoxic environments for extended periods, which may leave this species more vulnerable to demographic change.
An integrated approach that combines ecology, physiology and molecular tools will shed more light on the processes that control population patterns in this enigmatic group of fishes. This is especially important for understanding local adaptation, which may also vary considerably be-tween populations of the three species, especially within M. dorsalis. Based on our broad mtDNA study, genome-wide analyses are a future priority to determine levels of selection between populations of these fishes, which will help disentangle the relative contribution of gene flow and adaptation in shaping population structure in fishes with variable gene-flow patterns.