- Research article
- Open Access
Intraspecific genetic lineages of a marine mussel show behavioural divergence and spatial segregation over a tropical/subtropical biogeographic transition
BMC Evolutionary Biology volume 15, Article number: 100 (2015)
Intraspecific variability is seen as a central component of biodiversity. We investigated genetic differentiation, contemporary patterns of demographic connectivity and intraspecific variation of adaptive behavioural traits in two lineages of an intertidal mussel (Perna perna) across a tropical/subtropical biogeographic transition.
Microsatellite analyses revealed clear genetic differentiation between western (temperate) and eastern (subtropical/tropical) populations, confirming divergence previously detected with mitochondrial (COI) and nuclear (ITS) markers.
Gene flow between regions was predominantly east-to-west and was only moderate, with higher heterozygote deficiency where the two lineages co-occur. This can be explained by differential selection and/or oceanographic dynamics acting as a barrier to larval dispersal.
Common garden experiments showed that gaping (periodic closure and opening of the shell) and attachment to the substratum differed significantly between the two lineages. Western individuals gaped more and attached less strongly to the substratum than eastern ones.
These behavioural differences are consistent with the geographic and intertidal distributions of each lineage along sharp environmental clines, indicating their strong adaptive significance. We highlight the functional role of diversity below the species level in evolutionary trends and the need to understand this when predicting biodiversity responses to environmental change.
The early recognition of functionally significant levels of biodiversity within a species is a central objective of evolutionary biology and biodiversity conservation . This is particularly important at times when biodiversity loss is one of the planet´s major global problems [2,3]. Biodiversity can be investigated, and its importance evaluated, at different levels including species, their genes and functional traits [4,5].
Given the recurrent detection of distinct phylogeographic lineages, it is clearly necessary to consider the adaptive significance of intraspecific differences in order to understand the evolutionary potential of a species [6,7]. Complementary approaches based on neutral genetic markers, commonly employed to estimate recent demographic connectivity, and on methods that can reveal ecologically relevant adaptive traits across heterogeneous habitats are pivotal for the assessment and management of intraspecific diversity (e.g. [8,9]).
The identification of evolutionarily significant units uniquely based on neutral molecular markers, while overlooking intraspecific patterns of phenotypic divergence, can lead to underrepresentation of locally adapted populations and thus underestimate biodiversity [10,11]. Local environmental conditions resulting in spatially divergent selective regimes can maintain local adaptation across populations of the same species independently of varying and/or contrasting neutral genetic structure [12,13].
Species with wide distributions over variable environments and large effective population sizes often display marked intraspecific biodiversity with distinct genotypes/populations within a species, resulting from adaptation to local environmental conditions but capable of interbreeding with other ecotypes of the same species (e.g. [13,14]). Marine broadcast spawners with broad geographic ranges extending over abrupt environmental clines are ideal model organisms for studying the interplay between local adaptation and genetic structure due to neutral processes, because divergent selection produced by distinct environments is coupled with a high potential for gene flow among geographically separated subpopulations through larval transport.
The southern African coastline covers a wide range of climatic and oceanic conditions and can be divided into three major biogeographic regions: cool-temperate (from southern Namibia to the Cape of Good Hope; Figure 1A), warm-temperate (from the Cape of Good Hope eastward to East London) and sub-tropical (from East London north to Mozambique). The environmental discontinuities characteristic of these different parts of the coastline are mirrored by a degree of endemism among the different biogeographic provinces and phylogeographic breaks in species that occur in more than one province ( and references therein, ).
The brown mussel Perna perna is a habitat-forming species distributed widely around the world in warm waters. In southern Africa, it dominates intertidal shores in the sub-tropical and warm-temperate bioregions (from central Mozambique to the Cape of Good Hope, GH). It is absent from there to central Namibia due to the cold, upwelled waters of the Benguela system (a distributional gap of more than 1000 km; Figure 1A). It then extends northwards along the west coast of Africa to the southern Iberian Peninsula and into the Mediterranean Sea as far as the Gulf of Tunis ( and references therein).
Analyses of mitochondrial DNA (mtDNA) sequences indicate a sharp phylogeographic break on the southeast coast of South Africa (corresponding to the warm-temperate/subtropical transition) forming a western and an eastern lineage . The western lineage includes animals from both the Namibian coast and the south coast of South Africa (despite the wide gap across the Benguela upwelling system). The eastern lineage includes individuals from the southeast and east coasts of South Africa. The two lineages overlap in their distributions on the southeast coast over a distance of approximately 200 km. The presence of the two lineages was recently confirmed by nuclear (ITS) sequence data . Most importantly, mitochondrial (COI) and nuclear (ITS) markers clearly show a non-sister relationship for the two South African lineages . The most likely biogeographic scenario explaining the independent origins of the two South African lineages involves an Indo-Pacific origin for P. perna, dispersal into the Mediterranean and Atlantic through the Tethys seaway, and recent secondary contact after southward expansion of the western and eastern South African lineages . Manipulative experiments and the release of oceanographic drifters suggest that phylogeographic patterns of P. perna may be maintained by a combination of local conditions experienced in the different biogeographic regions, and by the isolating effect of regional oceanographic dynamics .
Here, we examine the hypotheses that contemporary demographic connectivity and functional traits of P. perna populations are influenced by the tropical/subtropical biogeographic transition. First, we adopted a population genetic approach using microsatellite markers, which are unaffected by natural selection, to investigate differentiation and migration patterns of P. perna between populations inhabiting different bioregions. The re-examination of genetic structure in South African P. perna with microsatellites and comparison with previously reported mitochondrial data is critical to assess potential mito-nuclear discordance and contemporary biogeographic conditions. Second, we used a common-garden design to compare fundamental behavioural traits known to affect survival along large (tropical/subtropical) and small (within-shore) scale environmental gradients. In particular, we assessed shell gaping and attachment strength, key attributes affecting tolerance of the two major stresses in intertidal habitats, desiccation and wave action. A common-garden design is necessary to distinguish between adaptive and plastic responses of population subjected to distinct selective pressures. Finally, we compared intertidal elevation (vertical distributional limits) of the two lineages by using in situ temperature profiles to compare intertidal limits of populations from the two bioregions.
Materials and methods
Mussel collection were performed under permit number RES2014/12 issued by the Department of Agriculture Forestry and Fisheries to the Department of Zoology and Entomology at Rhodes University.
DNA extraction, amplification and genotyping
Mussels (adults, >4cm in shell length) were collected between February 2010 and May 2013 at six locations in southern Africa (Figure 1A and Table 1) in the middle portion of the species’ vertical distribution. Mantle tissue (20-30 mg) was dissected from each individual, preserved in 96% ethanol and stored at -20°C. Total genomic DNA extraction was performed using a standard Proteinase K protocol adapted from . Molecular genotyping of a total of 269 samples was carried out using a set of 10 microsatellite loci in three multiplex reactions, with subsequent separation of the PCR products using an ABI PRISM 3130xl DNA analyser (Applied Biosystems) with GeneScan Liz 500 as size standard (Applied Biosystems) and following the procedure of Coelho et al. . Data were scored using PEAK-SCANNER v. 1.0 software (Applied Biosystems).
Scoring errors, large allele dropout and null alleles were checked employing the program MICROCHECKER  and frequency of null alleles estimated based on the algorithm presented in Brookfield . Additionally, tests for Hardy–Weinberg equilibrium were conducted using GENETIX 4.05 software  for all locus-location combinations by assessing the inbreeding estimator Wright’s FIS and estimating significant values after 10,000 permutations. Tests for linkage disequilibrium between all pairs of loci were performed according to the method of Black and Kraftsur  implemented in GENETIX.
For each location, allele frequencies, observed heterozygosity (HO) and expected heterozygosity (HE), plus FIS (significance based on 1,000 permutations and adjusted using q-value correction for multiple comparisons; ) were computed with GENETIX. Standardized allelic richness (Â) was estimated by resampling 1000 times to standardize to the smallest sample size (n = 30) to account for differences between sampled locations, using StandArich R package . Pairwise location estimates of genetic differentiation plus 95% bootstrapped confidence intervals (10,000 replicates) were estimated using Weir and Cockerham’s  FST estimator (θ) and Jost’s  DST. Calculations were executed in the Diversity R package .
To evaluate the extent of admixture between the two mitochondrial lineages using microsatellite data, we used STRUCTURE 2.2 . Clusters (K) varied from one to a maximum of seven, corresponding to the number of locations in the study data plus one. Twenty replicate runs [100,000 Markov Chain Monte Carlo (MCMC) iterations and 10,000 burn-in] were performed under the admixture and the correlated allele frequency model. The most probable value of K was inferred based on the ad hoc criteria L(K) proposed by Pritchard et al.  and delta K by Evanno et al. . Similar replicate runs were grouped based on the symmetric similarity coefficient of >0.9 using the FullSeach algorithm in CLUMPP 1.2  and visualized with bar_plotter.rb 1.1 implemented in Ruby.
To complement the results of STRUCTURE, we used discriminant analysis of principal components (DAPC, ), which can be performed when uncovering population admixture in the absence of assumptions. This method first transforms the data using principal components analysis, which ensures that the variables are not correlated and that the number of variables is smaller than the number of individuals. We ran two DAPC analyses in ADEGENET R package , following Jonker et al. : a location assignment and a genetic cluster assignment. We assigned each individual a priori to its location of origin and obtained for each one the probability of assignment to the correct sampling site. This procedure maximizes the among-location variation and minimizes the within-location variation . Using the find.clusters function (with 107 iterations), we ran successive K-means clustering of the individuals from K = 1 to K = 6, and identified the best supported number of clusters through comparison with the Bayesian Information Criterion (BIC) for the different values of K. We determined the number of clusters and assigned each individual to a genetic cluster without providing any a priori population assignment. In that way, the grouping factor is the genetic cluster and not the sampling location, offering an unbiased interpretation of population structure. To avoid unstable assignments of individuals to clusters, we retained 89 PCs (sample size divided by three), but used all discriminant functions, in a preliminary DAPC run. The results were then reiterated by the optim.a.score function with 100 simulations to determine the optimal number of PCs, and a final DAPC was subsequently carried out with the optimal number of PCs.
Estimating the directionality of gene flow
To compare different biogeographic hypotheses on migration rates of P. perna, we used the coalescent-based program Migrate-N MIGRATE version 3.1.3 [38,39]. This approach assumes the Wright-Fisher model, where locations have a constant effective size through time, the rate of mutation is constant, and locations exchange migrants with constant rates per generation, but those rates can vary among locations. We conducted the analyses on a random subsample of 30 individuals, with the dataset structured into three groups according to the genetic clusters recovered (i.e. all populations were assigned to one of the two clear clusters with the exception of PA, which fell indistinctly into each of the two, and so was considered separately): West (sites HM, PL), Centre (PA) and East (PE, DU, PO). We tested four variations of the three-group (West, Centre, East) migration model: all directional migration (Model 1: West ⟷ Centre ⟷ East, full model, nine parameters), strict west to east migration (Model 2: West → Centre → East, six parameters), strict east to west migration (Model 3: West ← Centre ← East, six parameters) and from the centre to the margins (Model 4: West ← Centre → East, five parameters). Testing the directionality of gene flow is justified because the dominant ocean current between the ocean basins, the Agulhas Current, runs westerly from the Indian towards the Atlantic Ocean (Model 3) and is thought to play a role in limiting marine dispersal in the opposite direction [40,41] though coastal counter-currents also occur  (Model 2). Models 1 and 4 consider indistinct migration among groups and migration directionality towards the transition area between the two clusters respectively. Initial values were calculated using F ST. Mutation rates were set to be constant among loci. The Migrate-N run parameters were calibrated on the full model for convergence of the Markov chain Monte Carlo sampling method. The prior distributions were uniform for mutation-scaled population size parameters theta (θ), that are four times the product of the effective population size and the mutation rate, and mutation-scaled migration rates M, that is, migration rate scaled by the mutation rate, over the range of θ = 0–100 and M = 0–100. These settings resulted in converged posterior distributions with a clear maximum for each estimate. The three replicate Bayesian runs consisted of one long chain with a total of 6 million states and genealogy changes after discarding the first 10,000 genealogies as “burn-in” for each locus. For all the analyses we used an adaptive heating scheme with 4 simultaneous chains using different acceptance ratios (temperature settings were 1.0; 1.5; 3.0; 1,000,000.0). Overall loci information was combined into a single estimate by Bézier approximation of the thermodynamic scores as described by Beerli and Palczewski  and we averaged the Bézier score over three different runs and used this as input to evaluate multiple models using Bayes factors .
Attachment strength and gaping behaviour
All mussels used for the assessment of behavioural traits were collected from locations known to support pure mtDNA lineages (). Specimens (shell length 3-4cm) were collected in February 2013 [at Brenton-on-Sea (BS) for the western lineage and at Amanzinoti (AM) for the eastern lineage; Figure 1A) and acclimated in seawater at 17°C for two weeks prior each experiment. A second trial was replicated in February 2015 with mussels from two additional locations [at Robberg (RO) for the western lineage and at Port St. Johns (PJ) for the eastern lineage; Figure 1A]. Common garden experiments involve the acclimation and comparison of individuals from distinct geographical locations under identical environmental conditions. Such experimental set-ups are widely used to disentangle the effects of genetic and environmental conditions on observed phenotype differences (e.g. ).
Mussels were divided in three separate tanks (n =5 of each lineage in each tank), where they were placed at least 10cm apart so that they maintained a solitary position (i.e. horizontal to the substratum) and attachment strengths were not influenced by nearby conspecifics. After three days, the attachment strength of each mussel was measured by drilling a hole into one of the valves and using a spring balance to record the force required to dislodge the mussel perpendicularly to the substratum .
Mussels (n = 15 for each lineage and each measuring interval) were exposed to air for 3h in a controlled environment chamber at 35°C and 60-80% humidity. For each one hour interval, the number of valve movements (gaping behaviour) were noted by visual observation, without recording the width of the gape .
Intertidal distributional limits
Temperature dataloggers (iButtons®, Maxim Integrated Products, Dallas Semiconductor, USA) were used to relate the effective shore height of P. perna at two open coast sites on the east coast [Amanzinoti (AM) and Balito (BA); Figure 1A] and two within the western range of the species (BS and RO). At each site, two transects were established, approximately 50 m apart and running perpendicular to the shoreline. Along each transect, two logger were deployed at the upper and two at the lower limit of P. perna distribution. Loggers were programmed to record data at 10 min intervals between February 1st to 14th of 2013 (western range) and March 17 to 30 of 2014 (eastern range), covering the two-week tidal cycle in both regions (i.e. a neap and a spring tide). Temperature logger profiles clearly reveal when loggers are first inundated by the returning tide (sudden, sharp temperature drops; ≥ 3°C over 10 min; ). This was used to calculate the effective shore height using the online forecasts service of SHOM database (http://www.shom.fr; settings: hauteur d’eau à une heure donnée). Because tide cycles are not in phase across large geographical scales, Durban and Mossel Bay tidal elevations were used for the east and south coast respectively. Where both loggers were recovered at a transect, the values were averaged and used to estimate populations’ intertidal limits at each transect.
Each trial was analysed separately. Attachment strength data were analysed using a 2-way ANOVA with Lineage (western or eastern) and Block (tank 1, 2 or 3) as fixed and random factors respectively. Gaping behaviour data were analysed using a 2-way ANOVA with Lineage (western, eastern) and Time (hour 1, 2 or 3) as fixed factors.
For each intertidal Limit (low, high), data were analysed under a nested design with Range (west, east) as a fixed factor, Site (1, 2) nested within Range and Transect nested within Range and Site.
Heteroscedasticity was tested using Levene’s test and post-hoc separation of significant means by SNK tests. Analyses were performed in SPSS (IBM Corp., USA).
Null alleles at high frequencies (>0.2) were indicated by MICROCHECKER for three loci (P11, P16, P27) for several localities (Additional file 1: Figure A1). For these three loci, departures from HWE were consistently reported across all localities. To avoid overestimation of FIS and FST estimators  by null alleles, these three loci were removed from all subsequent analyses, including descriptive analyses reported in Table 1. Particularly high null allele frequencies have often been reported in several invertebrate taxa, including corals and bivalves, due to high levels of DNA sequence variation . Specifically, nucleotide variations in microsatellite flanking regions can cause a high frequency of null-alleles, even within populations used for characterization of microsatellites .
All microsatellite loci were highly polymorphic among the six populations. The 269 individuals contained 226 alleles in these loci, with an average number of alleles per locus of 33, ranging from 9 (P29) to 98 (P20; Additional file 2: Figure A2). Allelic richness and gene diversity HE varied between 14.2-16.3 and 0.829-0.777 respectively and the highest allelic richness was recorded at DU (Table 1). Mean heterozygote deficiency was detected at all locations with the exception of PO (Table 1 and Additional file 3: Table A1 for estimates for each locus).
Genetic differentiation among populations was detected as indicated by pairwise FST and DST values (p < 0.05 for all comparisons with the exception of HM vs. PL and PE vs. PO; Table 2). While DST pairwise differentiation was lower within than among two groups (hereafter ‘western’ group: HM, PL; and ‘eastern’ group: PE, DU, PO), FST values did not clearly differentiate PA as belonging to either group.
Discriminant Analysis of Principal Components (DAPC), run with geographic location defined a priori, clearly differentiated the ‘western’ and ‘eastern’ groups on the scatter plot x-axis, confirming the presence of the genetic structure indicated by the F ST analysis (Figure 1B). Individuals from PA contained an admixture of alleles from both groups and indistinctly clustered with each of them. When the a posteriori genetic cluster assignment was run for K = 2 (best supported number of clusters), most individuals of the ‘western’ group were assigned to the same cluster (cluster 1; Figure 1C). The individuals of the ‘eastern’ group were also mostly assigned to one cluster (cluster 2), and in a higher rate than the populations of the first cluster. This is particularly true for PE, where all individuals were assigned to cluster 2. Individuals from PA were assigned to both clusters with a slightly higher number of individuals assigned to cluster 2.
This pattern of subdivisions was largely consistent with the results obtained with STRUCTURE, which revealed the most significant increase of ΔK at two clusters, (Figure 1D), thereafter ΔK remained unchanged. The L(K) plot rejected the option of K = 1. We concluded that K = 2 is the most likely number of genetic clusters (Additional file 4: Figure A3). While at most locations individuals displayed an admixture higher than 0.85 for either group or cluster, PA showed a split of 0.53 and 0.47 for the western and eastern groups respectively.
When comparing the four biogeographic hypotheses with MIGRATE-n, Model III, strictly east to west migration, was the only model supported (Table 3).
In each trial, lineage (as defined by mtDNA in ) was a highly significant factor influencing attachment strength (two-way ANOVA, square-root transformed, p < 0.01; Figure 2A), with the eastern lineage having a higher attachment strength than the western.
In each trail, individuals from the eastern lineage always gaped significantly more than those from the western lineage (two-way ANOVA, square-root transformed for the first trial only, p < 0.001; Figure 2B). Post-hoc SNK tests indicated that both lineages gaped significantly more at hour one than at hours two or three during the first trial (p < 0.001), however, during the second trail, there was a time × lineage interaction (p < 0.05) because at hour two eastern lineage individuals gaped more that at hour three.
P. perna in the western range reached a significantly higher intertidal elevation than on the east coast (three-way ANOVA, p < 0.001; n = 6; Figure 2C). However, mussels in the western and eastern ranges displayed the same lower vertical limit.
The microsatellite data revealed strong genetic division between the two geographically defined groups of populations, supporting previous findings based on mitochondrial and nuclear markers [18,19]. Populations from cluster one inhabit the warm-temperate, south-western coast of South Africa, while cluster two comprises populations distributed along the sub-tropical eastern shores. Most importantly, genetic admixture was only moderate, asymmetric and not restricted to the previously described 200 km contact zone, which lies mostly in the warm-temperate transition zone ).
The admixture detected with microsatellite markers is probably the result of secondary contact between ancient lineages; phylogeographic patterns based on mtDNA (COI) and nDNA (ITS) sequence data indicate an Indo-Pacific origin of the species and a dispersal into the Mediterranean and Atlantic through the Tethys seaway. After the closing of the seaway, the two non-sister lineages diverged (310 kyr) and came into secondary contact after their respective eastern, or western, southward expansion . However, genetic differentiation between lineages can be primarily attributed to the magnitude of the differences in allele frequencies rather than to unique alleles, as would be expected in ancient isolated lineages.
The detected gene exchange was primarily unidirectional, from east to west. Additionally, the sampled location at the geographic border of the two lineages (PA) shows the strongest gene intermediacy. Gene exchange and the maintenance of this pattern could be attributed to several non-mutually exclusive effects. Eastward gene flow could be impeded by dispersal barriers (e.g. currents), or by selection (i.e. unfavourable environmental conditions) causing outbreeding depression if lineages are locally adapted. Alternatively, poor mating success between lineages could also play a role, either due to divergent reproductive ecology, or to some level of genetic incompatibility.
Barriers to dispersal
Mussel larval dispersal is largely influenced by local hydrodynamics. The main oceanographic influence on the east and south coasts of South Africa is the Agulhas Current. This is a major ocean current that flows towards the southeast, bringing warm, oligotrophic water from the Mozambique channel . The inshore thermal front of this current varies geographically and can bring its warm waters very close inshore [53,54]. From near PA southwards the current is deflected away from the coast. The trajectories of nearshore drogues released off the south coast revealed no mixing of southern waters into the eastern water mass, while drifters released off the east coast moved southward, caught up in the Agulhas Current . Although there are wind-driven surface currents that could allow limited eastward larval dispersal , water exchange is predominantly unidirectional, promoting westerly larval dispersal of the eastern lineage and explaining the asymmetrical gene flow described here. Previous mtDNA data on several species with planktonic larvae (including P. perna; ) support our results obtained with microsatellite data, by showing asymmetric east-to-west migration along the east coast (but see ).
Interestingly, one of our most eastern locations (DU) displayed a fair amount of genetic admixture. DU is the largest and busiest shipping terminal on the African Continent, and globally one of the seven key epicentres for inter-regional exchange of species . Shipping activities are known vectors of larval dispersal [58,59]. The fact that DU is the only location on the east coast with these characteristics suggests that the introduction of non-indigenous genes from the western range detected only here is a consequence of human-mediated larval transport.
Departures from random mating
Out of the six populations sampled, five displayed significant heterozygote deficiency. Null alleles are unlikely to be the cause of high Fis values in this case, since microsatellite loci with a high frequency of null alleles were excluded from the analyses, while the remaining loci compiled in the final dataset displayed similar homozygote excess simultaneously across all loci. Inbreeding is also unlikely in a broadcast spawner with long a larval pelagic phase (assumed to be similar to that of Mytilus spp., [60,61]) as recruitment is expected to produce spatially random kinship.
One possible explanation for heterozygote deficiency is a Wahlund effect due to the presence of genetically distinct subpopulations (which may be in HW equilibrium) in each population sampled. In this context, it is interesting that PA, which is at the contact zone between the two lineages, displayed the lowest mean observed heterozygosity, although this was not true across all loci (see Additional file 3: Table A1). Such effects can be expected when pooling discrete subpopulations with different allele frequencies that do not interbreed as a single randomly mating unit. In contrast, non-significant levels of heterozygote deficiency were detected in PO, which is the farthest population from the contact zone and located on the east side of the coast.
Reproductive isolation of the lineages, leading to heterozygote deficiency, could result from differences in the timing and frequency of spawning events [46,62,63]. It is not known when each lineage spawns nor what triggers the event, but lack of synchrony in spawning by co-existing mussel species has been reported, and analogous intraspecific differences may occur between lineages of P. perna .
Behaviour can radically moderate an organism’s experience of the environment and so dictate physiological reactions, and consequently, tolerance limits to the stress imposed by those environmental conditions. We show that the two P. perna mtDNA lineages exhibit clear differences in behavioural traits, indicating an intraspecific behavioural divergence that potentially influences vertical zonation and habitat segregation over large and small spatial scales. By acclimating and comparing individuals from different locations in a common laboratory environment (common garden) we could differentiate between adaptive and plastic responses to distinct environmental conditions. The maintenance of behavioural phenotypic divergence in experimental populations indicated that differences are due to underlying genetic divergence. Mussel gaping, the intermittent opening and closure of the valves, allows the maintenance of aerobic respiration by sustaining an oxygen gradient across the mantle wall and the gills during emersion. Respiration rates of organisms usually increase with temperature and several studies have shown that oxygen consumption by mussels increases exponentially with temperature (e.g. ). In P. perna, gaping significantly increased with temperature, supporting the important role of this behaviour in aerobic respiration  and the correlation between gaping and heat stress.
Although gaping has advantageous respiratory effects, it results in water loss caused by both evaporation and incidental discharge of water during valve closure, thus increasing the risk of desiccation . Although valve opening suggests possible evaporative cooling, several studies have excluded any link between this behaviour and the body temperature of isolated individuals [66,67] (but see  for emergent cooling effects in mussel aggregations).
The adaptive significance of gaping behaviour helps explain the vertical intertidal zonation of the western lineage of P. perna and the coexisting invasive Mytilus galloprovincialis, a non-gaping mussel species [47,68] highlighting how even very simple behaviour can influence the outcome of interactions between coexisting species.
Our results and previous studies [47,68] suggest a determinant effect of gaping on the small and large scale spatial distributions of P. perna lineages. More frequent gaping in the eastern lineage would promote efficient respiration under high, subtropical air temperatures, but more gaping would simultaneously increase the risk of desiccation, restricting the eastern lineage to lower on the shore than western conspecifics ([as has been shown for a comparison of P. perna with M. galloprovincialis; ).
In addition to differences in gaping rates, the eastern lineage showed much stronger attachment to the substratum. Firm byssal attachment is a prerequisite for survival in the intertidal (e.g. [69,70]) but this difference is unlikely to explain large scale spatial segregation of the lineages as the temperate and subtropical regions are not known to differ significantly in terms of wave action.
Temperature offers a more likely explanation as it has direct or indirect effects on attachment. Higher temperatures promote byssal thread production while weakening attachment by accelerating thread deterioration and negatively affecting the curing and moulding of byssal threads [71,72]. Additionally, collagenases are common in marine bacteria  and their activity is enhanced by temperature , potentially contributing to the weakening of the byssus. Thus, the greater attachment strength of mussels from the east coast could be an adaptive response to the thermally more challenging subtropical environment.
This study shows how distinct genetic lineages, partially segregated along a sharp geographic environmental gradient display adaptive responses to different selective pressures. P. perna inhabits shores that are experiencing rapid and geographically uneven environmental changes [75,76]. Our results clearly indicate that, rather than the species exhibiting a common overall response, the two lineages are likely to be affected by, and react quite differently to, climate change. Our findings highlight how critical it is to understand intraspecific diversity if we wish to have a grasp of the real magnitude of biodiversity and how within-species diversity will shape ecosystem resilience/resistance in a changing environment.
Availability of supporting data
The data set supporting the results of this article is included within the article (Additional file 5: Table A2).
Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM. Beyond predictions: biodiversity conservation in a changing climate. Science. 2011;332(6025):53–8.
Stokstad E. Despite progress, biodiversity declines. Science. 2010;329(5997):1272–3.
Cardinale BJ, Duffy JE, Gonzalez A, Hooper DU, Perrings C, Venail P, et al. Biodiversity loss and its impact on humanity. Nature. 2012;486(7401):59–67.
Reusch TBH, Ehlers A, Hämmerli A, Worm B. Ecosystem recovery after climatic extremes enhanced by genotypic diversity. Proc Natl Acad Sci U S A. 2005;102(8):2826–31.
Cardinale BJ, Wright JP, Cadotte MW, Carroll IT, Hector A, Srivastava DS, et al. Impacts of plant diversity on biomass production increase through time because of species complementarity. Proc Natl Acad Sci. 2007;104(46):18123–8.
Violle C, Enquist BJ, McGill BJ, Jiang L, Albert CH, Hulshof C, et al. The return of the variance: intraspecific variability in community ecology. Trends Ecol Evol. 2012;27(4):244–52.
Nicastro KR, Zardi GI, Teixeira S, Neiva J, Serrao EA, Pearson GA. Shift happens: trailing edge contraction associated with recent warming trends threatens a distinct genetic lineage in the marine macroalga Fucus vesiculosus. BMC Biol. 2013;11(1):6.
Teske P, Papadopoulos I, Newman B, Dworschak P, McQuaid C, Barker N. Oceanic dispersal barriers, adaptation and larval retention: an interdisciplinary assessment of potential factors maintaining a phylogeographic break between sister lineages of an African prawn. BMC Evol Biol. 2008;8(1):341.
Ballentine B, Greenberg R. Common garden experiment reveals genetic control of phenotypic divergence between Swamp Sparrow subspecies that lack divergence in neutral genotypes. PLoS ONE. 2010;5(4):e10229.
Zink RM. The role of subspecies in obscuring avian biological diversity and misleading conservation policy. Proc R Soc Lond Ser B Biol Sci. 2004;271(1539):561–4.
Phillimore AB, Owens IPF. Are subspecies useful in evolutionary and conservation biology? Proc R Soc B Biol Sci. 2006;273(1590):1049–53.
Crispo E. Modifying effects of phenotypic plasticity on interactions among natural selection, adaptation and gene flow. J Evol Biol. 2008;21(6):1460–9.
Zardi GI, Nicastro KR, Ferreira Costa J, Serrão EA, Pearson GA. Broad scale agreement between intertidal habitats and adaptive traits on a basis of contrasting population genetic structure. Estuar Coast Shelf Sci. 2013;131:140–8.
Sambatti JBM, Rice KJ. Local adaptation, patterns of selection, and gene flow in the californian serpentine sunflower (Helianthus exilis). Evolution. 2006;60(4):696–710.
Awad AA, Griffiths CL, Turpie JK. Distribution of South African marine benthic invertebrates applied to the selection of priority conservation areas. Divers Distrib. 2002;8(3):129–45.
Teske PR, Von Der Heyden S, McQuaid CD, Barker NP. A review of marine phylogeography in southern Africa. S Afr J Sci. 2011;107(5-6):43–53.
Lourenço C, Nicastro KR, Serrão EA, Zardi GI. First record of the brown mussel (Perna perna) from the European Atlantic coast. Marine Biodiversity Record. 2012;5:e39.
Zardi G, McQuaid C, Teske P, Barker N. Unexpected genetic structure of mussel populations in South Africa: indigenous Perna perna and invasive Mytilus galloprovincialis. Mar Ecol Prog Ser. 2007;337:135–44.
Cunha RL, Nicastro KR, Costa J, McQuaid CD, Serrão EA, Zardi GI. Wider sampling reveals a non-sister relationship for geographically contiguous lineages of a marine mussel. Ecology Evol. 2014;4(11):2070–81.
Zardi G, Nicastro K, McQuaid C, Hancke L, Helmuth B. The combination of selection and dispersal helps explain genetic structure in intertidal mussels. Oecologia. 2011;165(4):947–58.
Sambrook J, Fritseh EF, Maniatis T. Molecular Cloning. A Laboratory Manual. 2nd ed. Cold Spring Harbor, NY: Cold Spring Harbor Laboratory Press; 1989.
Coelho N, Zardi GI, Pearson GA, Serrao EA, Nicastro KR. Characterization of ten highly polymorphic microsatellite loci for the intertidal mussel Perna perna, and cross species amplification within the genus. BMC Res Notes. 2012;5(1):558.
Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. Micro-checker: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4(3):535–8.
Brookfield JFY. A simple new method for estimating null allele frequency from heterozygote deficiency. Mol Ecol. 1996;5(3):453–5.
Belkhir, K., et al. GENETIX 4.05, logiciel sous Windows TM pour la génétique des populations. Laboratoire génome, populations, interactions, CNRS UMR 5000 1996;1996-2004.
Black WC, Krafsur ES. A FORTRAN program for the calculation and analysis of two-locus linkage disequilibrium coefficients. TAG Theor Applied Genetics Theoretische und angewandte Genetik. 1985;70(5):491–6.
Storey JD. A direct approach to false discovery rates. J R Stat Soc Ser B. 2002;64:479–98.
Alberto F, Arnaud-Haond S, Duarte CM, Serrão EA. Genetic diversity of a clonal angiosperm near its range limit: the case of Cymodocea nodosa at the Canary Islands. Marine Ecology-Progress Series. 2006;309:117–29.
Weir BS, Cockerham CC. Estimation of gene flow from F-statistics. Evolution. 1984;47:855–63.
Jost LOU. GST and its relatives do not measure differentiation. Mol Ecol. 2008;17(18):4015–26.
Keenan K, McGinnity P, Cross TF, Crozier WW, Prodöhl PA. diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol Evol. 2013;4(8):782–8.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software structure: a simulation study. Mol Ecol. 2005;14(8):2611–20.
Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23(14):1801–6.
Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010;11(1):94.
Jombart T, Devillard S, Dufour AB, Pontier D. Revealing cryptic spatial patterns in genetic variability by a new multivariate method. Heredity. 2008;101(1):92–103.
Jonker R, Kraus R, Zhang Q, Hooft P, Larsson K, Jeugd H, et al. Genetic consequences of breaking migratory traditions in barnacle geese Branta leucopsis. Mol Ecol. 2013;22(23):5835–47.
Beerli P, Felsenstein J. Maximum likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genetics. 1999;152:763–73.
Beerli P, Felsenstein J. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proc Natl Acad Sci. 2001;98(8):4563–8.
Ivanova E. The global thermohaline paleocirculation. Springer Science & Business Media, 2009.
Lutjeharms J. The agulhas current. Afr J Mar Sci. 2006;28:729–32.
Jackson JM, Rainville L, Roberts MJ, McQuaid CD, Lutjeharms JRE. Mesoscale bio-physical interactions between the Agulhas Current and the Agulhas Bank, South Africa. Cont Shelf Res. 2012;49:10–24.
Beerli P, Palczewski M. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics. 2010;185(1):313–26.
Bloomquist E, Lemey P, Suchard M. Three roads diverged? Routes to phylogeographic inference. Trends Ecol Evol. 2010;25(11):626–32.
Pearson GA, Lago-Leston A, Mota C. Frayed at the edges: selective pressure and adaptive response to abiotic stressors are mismatched in low diversity edge populations. J Ecol. 2009;97(3):450–62.
Zardi GI, McQuaid CD, Nicastro KR. Balancing survival and reproduction: seasonality of attachment strength and reproductive output in indigenous (Perna perna) and invasive (Mytilus galloprovincialis) mussels. Mar Ecol Prog Ser. 2007;334:155–67.
Nicastro KR, Zardi GI, McQuaid CD, Stephens L, Radloff S, Blatch GL. The role of gaping behaviour in habitat partitioning between coexisting intertidal mussels. BMC Ecol. 2010;10(1):17.
Harley CDG, Helmuth BST. Local- and regional- scale effects of wave exposure, thermal stress, and absolute versus effective shore level on patterns of intertidal zonation. Limnol Oceanogr. 2003;48(4):1498–508.
Chapuis M-P, Estoup A. Microsatellite null Alleles and estimation of population differentiation. Mol Biol Evol. 2007;24(3):621–31.
Selkoe KA, Toonen RJ. Microsatellites for ecologists: a practical guide to using and evaluating microsatellite markers. Ecol Lett. 2006;9(5):615–29.
Hedgecock D, Li G, Hubert S, Bucklin K, Ribes V. Widespread null alleles and poor cross-species amplification of microsatellite DNA loci cloned from the Pacific oyster, Crassostrea gigas. J Shellfish Res. 2004;23(2):379–86.
Lutjeharms, Johann RE. The Agulhas Current. 2006:729-730.
Goschen WS, Schumann EH. An Agulhas Current intrusion into Algoa Bay during August 1988. S Afr J Mar Sci. 1994;14(1):47–57.
Schumann EH, Churchill JRS, Zaayman HJ. Oceanic variability in the western sector of Algoa Bay, South Africa. Afr J Mar Sci. 2005;27(1):65–80.
McQuaid CD, Phillips TE. Limited wind-driven dispersal of intertidal mussel larvae: in situ evidence from the plankton and the spread of the invasive species Mytilus galloprovincialis in South Africa. Mar Ecol Prog Ser. 2000;201:211–20.
Teske P, Papadopoulos I, Zardi G, McQuaid C, Edkins M, Griffiths C, et al. Implications of life history for genetic structure and migration rates of southern African coastal invertebrates: planktonic, abbreviated and direct development. Mar Biol. 2007;152(3):697–711.
Drake JM, Lodge DM. Global hot spots of biological invasions: evaluating options for ballast–water management. Proceedings of the Royal Society of London Series B: Biological Sciences 2004, 271(1539):575–580.
Carlton JT. 13 The scale and ecological consequences of biological invasions in the World's oceans. Invasive species and biodiversity management 2001, 24:195.
Flagella MM, Soria A, Buia MC. Shipping traffic and introduction of non-indigenous organisms: Study case in two Italian harbours. Ocean & Coastal Management 2006, 49(12):947–960.
Bayne BL. Reproduction in bivalve molluscs under environmental stress. Columbia: University of South Carolina Press; 1975.
de Vooys CGN. Numbers of larvae and primary plantigrades of the mussel Mytilus edulis in the western Dutch Wadden Sea. J Sea Res. 1999;41(3):189–201.
Nicastro KR, Zardi GI, McQuaid CD. Differential reproductive investment, attachment strength and mortality of invasive and indigenous mussels across heterogeneous environments. Biol Invasions. 2010;12(7):2165–77.
Lasiak T. The reproductive cycles of the intertidal bivalves Crassostrea cucullata (Born, 1778) and Perna perna (Linnaeus, 1758) from the Transkei Coast, southern Africa. Veliger. 1986;29:226–30.
Widdows J, Bayne BL, Livingstone DR, Newell RIE, Donkin P. Physiological and biochemical responses of bivalve molluscs to exposure in air. Comp Biochem Physiol. 1979;62:301–8.
Jansen JM, Hummel H, Bonga SW. The respiratory capacity of marine mussels (Mytilus galloprovincialis) in relation to the high temperature threshold. Comp Biochem Physiol A Mol Integr Physiol. 2009;153(4):399–402.
Fitzhenry T, Halpin PM, Helmuth B. Testing the effects of wave exposure, site, and behavior on intertidal mussel body temperatures: applications and limits of temperature logger design. Mar Biol. 2004;145(2):339–49.
Lent CM. Adaptations of the ribbed mussel, Modiolus demissus (Dillvvyn), to the intertidal habitat. Am Zoology. 1969;9:283–92.
Nicastro KR, Zardi GI, McQuaid CD, Pearson GA, Serrão EA. Love thy neighbour: group properties of gaping behaviour in mussel aggregations. PLoS ONE. 2012;7(10):e47382.
Zardi GI, Nicastro KR, McQuaid CD, Rius M, Porri F. Hydrodynamic stress and habitat partitioning between indigenous (Perna perna) and invasive (Mytilus galloprovincialis) mussels: constraints of an evolutionary strategy. Mar Biol. 2006;150(1):79–88.
Babarro J,MF, Abad MJ. Co-existence of two mytilid species in a heterogeneous environment: mortality, growth and strength of shell and byssus attachment. Mar Ecol Prog Ser. 2013;476:115–28.
Moeser GM, Carrington E. Seasonal variation in mussel byssal thread mechanics. J Exp Biol. 2006;209(10):1996–2003.
Garner YL, Litvaitis MK. Effects of wave exposure, temperature and epibiont fouling on byssal thread production and growth in the blue mussel, Mytilus edulis, in the Gulf of Maine. J Exp Mar Biol Ecol. 2013;446:52–6.
Merkel JR, Dreisbach JH, Ziegler HB. Collagenolytic activity of some marine bacteria. Appl Microbiol. 1975;29(2):145–51.
Bălan D, Israel-Roming F, Cornea P, Gherghina E, Luţă G, Matei F, et al. Novel fungal collagenase from Aspergillus oryzae. Sci Bull Ser F Biotechnol. 2013;17:160–3.
Lima FP, Wethey DS. Three decades of high-resolution coastal sea surface temperatures reveal more than warming. Nat Commun. 2012;3:704.
Teske PR, Zardi GI, McQuaid CD, Nicastro K. Two sides of the same coin: extinctions and originations across the Atlantic/Indian Ocean boundary as a consequence of the same climate oscillation. Frontiers Biogeography. 2013;5:48–59.
Emanuel B, Bustamante R, Branch G, Eekhout S, Odendaal F. A zoogeographic and functional approach to the selection of marine reserves on the west coast of South Africa. S Afr J Mar Sci. 1992;12(1):341–54.
This research was funded by project PTDC/BIA-BEC/103916/2008 from the Fundação para a Ciência e a Tecnologia (FCT; http://www.fct.pt/index.phtml.pt; to GIZ) and supported by an award from the South Africa Research Chairs Initiative (SARChI) of the Department of Science and Technology (to CDM). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We are grateful to R. Jacinto for technical support and P. Teske for kindly supplying specimens from Ponta Do Ouro.
There are no financial or non-financial competing interests.
GIZ, KRN: conceived and designed the experiment. JC, GIZ, KRN, CDM: performed the experiments. KRN, RC, JC: analyzed the data. GIZ, EAS, GAP: contributed reagents/materials/analysis tools. GIZ, KRN, CDM, GAP, EAS, RC: contributed to the writing of the manuscript. All authors read and approved the final manuscript.
Gerardo I Zardi and Katy R Nicastro contributed equally to this work.
Frequencies of null alleles per locus per population. Estimates based on the algorithm presented in Brookfield  by MICROCHECKER.
Allele frequencies of each location. Each box represents one locus, according to the following order: P01, P05, P11, P27, P29, P02, P08, P16, P20, P26.
Observed heterozygosity calculated for each locations and locus. Estimates made in GENETIX program .
Genotyping Table. Scored alleles per location and per locus.