Skip to content


BMC Evolutionary Biology

Open Access

Intraspecific genetic lineages of a marine mussel show behavioural divergence and spatial segregation over a tropical/subtropical biogeographic transition

  • Gerardo I Zardi1,
  • Katy R Nicastro2Email author,
  • Christopher D McQuaid1,
  • Rita Castilho2,
  • Joana Costa2, 3,
  • Ester A Serrão2 and
  • Gareth A Pearson2
Contributed equally
BMC Evolutionary Biology201515:100

Received: 17 December 2014

Accepted: 29 April 2015

Published: 31 May 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.


Genetic diversityMarine connectivity Perna perna


The early recognition of functionally significant levels of biodiversity within a species is a central objective of evolutionary biology and biodiversity conservation [1]. 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 ([15] and references therein, [16]).
Figure 1

Study area and genetic results. A) Map of the study area and sampling sites. Black dots: localities either mentioned in the text or used for ecologic data; blue dots: localities used for genetic data abbreviated as in Table 1. Continuous, black line: distribution of the western lineage (according to [18]); dashed, black line: distribution of the eastern lineage (according to [18]). The three biogeographic regions of southern Africa (according to [77]) are delimited by grey dotted lines. B) Results of the Bayesian approach using the Discriminant Analysis of Principal Components (DAPC) scatter plot of individuals with a priori-defined locations, to investigate genetic structure. Labels were placed at the centre of dispersion of each group, which are delineated by inertia ellipses. C) DAPC a posteriori genetic cluster assignment for K = 2. The size of each square is proportional to the percentage of individuals assigned to each cluster at each location. D) Stacked bar graph of assignment probabilities per individual within P. perna locations across the southern tip of South Africa. Assignment proportions for each individual were obtained from Structure, based on genotypes from seven nuclear microsatellite loci, for K = 2.

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 ([17] 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 [18]. 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 [19]. Most importantly, mitochondrial (COI) and nuclear (ITS) markers clearly show a non-sister relationship for the two South African lineages [19]. 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 [19]. 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 [20].

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

Ethics statement

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.

Genetic data

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 [21]. 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. [22]. Data were scored using PEAK-SCANNER v. 1.0 software (Applied Biosystems).
Table 1

Genetic diversity of each population








Â (±SD)



−34.41, 19.26





14.36 (±0.55)

Plettenberg Bay


−34.00, 23.46





14.2 (±0.16)

Port Alfred


−33.48, 27.15





14.74 (±0.29)

Port Edward


−31.05, 30.23





15.7 (±0.45)



−29.55, 31.22





16.3 (±0.17)

Ponta do Ouro


−26.84, 32.89





15.6 (±0)

Locations are depicted in the map of Figure 1A and are ordered from west to east. Decimal coordinates (latitude and longitude, respectively); N: sample sizes; HE and HO: expected and observed heterozygosity; FIS: inbreeding coefficient with significant values after correction for multiple test indicated with *(p < 0.05); Â (±SD): allelic richness (standardized for N = 30 by resampling) ± standard deviation.

Location-based analyses

Scoring errors, large allele dropout and null alleles were checked employing the program MICROCHECKER [23] and frequency of null alleles estimated based on the algorithm presented in Brookfield [24]. Additionally, tests for Hardy–Weinberg equilibrium were conducted using GENETIX 4.05 software [25] 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 [26] 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; [27]) 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 [28]. Pairwise location estimates of genetic differentiation plus 95% bootstrapped confidence intervals (10,000 replicates) were estimated using Weir and Cockerham’s [29] FST estimator (θ) and Jost’s [30] DST. Calculations were executed in the Diversity R package [31].

Estimating admixture

To evaluate the extent of admixture between the two mitochondrial lineages using microsatellite data, we used STRUCTURE 2.2 [32]. 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. [32] and delta K by Evanno et al. [33]. Similar replicate runs were grouped based on the symmetric similarity coefficient of >0.9 using the FullSeach algorithm in CLUMPP 1.2 [34] 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, [35]), 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 [36], following Jonker et al. [37]: 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 [35]. 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 [42] (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 [43] and we averaged the Bézier score over three different runs and used this as input to evaluate multiple models using Bayes factors [44].

Ecological data

Attachment strength and gaping behaviour

All mussels used for the assessment of behavioural traits were collected from locations known to support pure mtDNA lineages ([18]). 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. [45]).

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 [46].

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 [47].

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; [48]). This was used to calculate the effective shore height using the online forecasts service of SHOM database (; 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.

Data analyses

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).


Genetic data

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 [49] 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 [50]. Specifically, nucleotide variations in microsatellite flanking regions can cause a high frequency of null-alleles, even within populations used for characterization of microsatellites [51].

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.
Table 2

Genetic differentiation between pairs of populations















(-0.018, 0.032)

(0.024, 0.113)

(0.216, 0.377)

(0.142, 0.262)

(0.205, 0.381)








(-0.007, 0.008)


(0.029, 0.132)

(0.237, 0.404)

(0.142, 0.271)

(0.228, 0.413)








(0.03, 0.068)

(0.04, 0.081)


(0.142, 0.281)

(0.064, 0.187)

(0.157, 0.355)








(0.06, 0.1)

(0.071, 0.108)

(0.033, 0.065)


(0.021, 0.127)

(-0.018, 0.08)








(0.057, 0.09)

(0.067, 0.098)

(0.016, 0.048)

(0.01, 0.033)


(0.01, 0.118)









(0.067, 0.114)

(0.048, 0.105)

(-0.007, 0.032)

(0.015, 0.053)

Codes correspond to locations in Figure 1A and are ordered from west to east. Pairwise genetic differentiation (FST, below the diagonal, DST above the diagonal). 95% confidence intervals are shown in brackets and confidence intervals not enclosing 0 are represented in bold.

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).
Table 3

Bayes factors model comparison of migration models for P. perna


Model parameters




Set 1

Set 2

Set 3

Model 1

W↔ C ↔ E and W ↔ E







Model 2

W → C → E and W → E







Model 3

W ← C ← E and W ← E







Model 4

W ← C → E







Model parameters code as follows: the first set of three signs indicate estimation of theta for West, migration from Centre to West, migration from East to West; the second set indicate migration from West to Centre, estimation of theta for Centre, migration from East to Centre; the third set indicates migration from West to East, migration from Centre to East and estimation of theta for East. An asterisk indicates that that particular parameter was estimated by the model, and a 0 indicates that no migration was allowed.

Ecological data

In each trial, lineage (as defined by mtDNA in [18]) 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.
Figure 2

Ecological results. A) Gaping behaviour of the eastern and western lineages, each trial is plotted separately, B) Strength of attachment (blocks pooled), each trial plotted separately, C) High and low intertidal limits of P. perna western and eastern lineages. Dots depict mild outliers (i.e. outliers outside the inner fences) and stars are extreme outliers (i.e. outliers outside the outer fences).

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 [18]).

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 [19]. 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 [52]. 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 [20]. Although there are wind-driven surface currents that could allow limited eastward larval dispersal [55], 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; [56]) support our results obtained with microsatellite data, by showing asymmetric east-to-west migration along the east coast (but see [8]).

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 [57]. 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 [46].


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. [64]). In P. perna, gaping significantly increased with temperature, supporting the important role of this behaviour in aerobic respiration [65] 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 [47]. 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 [68] 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; [47]).

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 [73] and their activity is enhanced by temperature [74], 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).




This research was funded by project PTDC/BIA-BEC/103916/2008 from the Fundação para a Ciência e a Tecnologia (FCT;; 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.

Authors’ Affiliations

Centre of Marine Sciences – CCMAR, Campus de Gambelas, Universidade do Algarve, Faro, Portugal
Present address: Algal Biology Lab, School of Botany, The University of Melbourne, Parkville, Australia


  1. Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM. Beyond predictions: biodiversity conservation in a changing climate. Science. 2011;332(6025):53–8.View ArticlePubMedGoogle Scholar
  2. Stokstad E. Despite progress, biodiversity declines. Science. 2010;329(5997):1272–3.View ArticlePubMedGoogle Scholar
  4. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  5. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  6. 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.View ArticlePubMedGoogle Scholar
  7. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  8. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  9. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  10. 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.View ArticleGoogle Scholar
  11. Phillimore AB, Owens IPF. Are subspecies useful in evolutionary and conservation biology? Proc R Soc B Biol Sci. 2006;273(1590):1049–53.View ArticleGoogle Scholar
  12. Crispo E. Modifying effects of phenotypic plasticity on interactions among natural selection, adaptation and gene flow. J Evol Biol. 2008;21(6):1460–9.View ArticlePubMedGoogle Scholar
  13. 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.View ArticleGoogle Scholar
  14. 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.View ArticlePubMedGoogle Scholar
  15. 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.View ArticleGoogle Scholar
  16. 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.Google Scholar
  17. 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.View ArticleGoogle Scholar
  18. 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.View ArticleGoogle Scholar
  19. 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.Google Scholar
  20. 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.View ArticlePubMedGoogle Scholar
  21. Sambrook J, Fritseh EF, Maniatis T. Molecular Cloning. A Laboratory Manual. 2nd ed. Cold Spring Harbor, NY: Cold Spring Harbor Laboratory Press; 1989.Google Scholar
  22. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  23. 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.View ArticleGoogle Scholar
  24. Brookfield JFY. A simple new method for estimating null allele frequency from heterozygote deficiency. Mol Ecol. 1996;5(3):453–5.View ArticlePubMedGoogle Scholar
  25. 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.Google Scholar
  26. 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.View ArticleGoogle Scholar
  27. Storey JD. A direct approach to false discovery rates. J R Stat Soc Ser B. 2002;64:479–98.View ArticleGoogle Scholar
  28. 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.View ArticleGoogle Scholar
  29. Weir BS, Cockerham CC. Estimation of gene flow from F-statistics. Evolution. 1984;47:855–63.Google Scholar
  30. Jost LOU. GST and its relatives do not measure differentiation. Mol Ecol. 2008;17(18):4015–26.View ArticlePubMedGoogle Scholar
  31. 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.View ArticleGoogle Scholar
  34. 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.View ArticlePubMedGoogle Scholar
  35. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  36. 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.View ArticlePubMedGoogle Scholar
  37. 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.View ArticlePubMedGoogle Scholar
  38. 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.PubMed CentralPubMedGoogle Scholar
  39. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  40. Ivanova E. The global thermohaline paleocirculation. Springer Science & Business Media, 2009.Google Scholar
  41. Lutjeharms J. The agulhas current. Afr J Mar Sci. 2006;28:729–32.View ArticleGoogle Scholar
  42. 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.View ArticleGoogle Scholar
  43. Beerli P, Palczewski M. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics. 2010;185(1):313–26.View ArticlePubMed CentralPubMedGoogle Scholar
  44. Bloomquist E, Lemey P, Suchard M. Three roads diverged? Routes to phylogeographic inference. Trends Ecol Evol. 2010;25(11):626–32.View ArticlePubMed CentralPubMedGoogle Scholar
  45. 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.View ArticleGoogle Scholar
  46. 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.View ArticleGoogle Scholar
  47. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  48. 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.View ArticleGoogle Scholar
  49. Chapuis M-P, Estoup A. Microsatellite null Alleles and estimation of population differentiation. Mol Biol Evol. 2007;24(3):621–31.View ArticlePubMedGoogle Scholar
  50. Selkoe KA, Toonen RJ. Microsatellites for ecologists: a practical guide to using and evaluating microsatellite markers. Ecol Lett. 2006;9(5):615–29.View ArticlePubMedGoogle Scholar
  51. 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.Google Scholar
  52. Lutjeharms, Johann RE. The Agulhas Current. 2006:729-730.Google Scholar
  53. Goschen WS, Schumann EH. An Agulhas Current intrusion into Algoa Bay during August 1988. S Afr J Mar Sci. 1994;14(1):47–57.View ArticleGoogle Scholar
  54. 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.View ArticleGoogle Scholar
  55. 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.View ArticleGoogle Scholar
  56. 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.View ArticleGoogle Scholar
  57. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  58. Carlton JT. 13 The scale and ecological consequences of biological invasions in the World's oceans. Invasive species and biodiversity management 2001, 24:195.Google Scholar
  59. 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.View ArticleGoogle Scholar
  60. Bayne BL. Reproduction in bivalve molluscs under environmental stress. Columbia: University of South Carolina Press; 1975.Google Scholar
  61. 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.View ArticleGoogle Scholar
  62. 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.View ArticleGoogle Scholar
  63. 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.Google Scholar
  64. 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.View ArticleGoogle Scholar
  65. 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.View ArticlePubMedGoogle Scholar
  66. 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.View ArticleGoogle Scholar
  67. Lent CM. Adaptations of the ribbed mussel, Modiolus demissus (Dillvvyn), to the intertidal habitat. Am Zoology. 1969;9:283–92.Google Scholar
  68. 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.View ArticlePubMed CentralPubMedGoogle Scholar
  69. 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.View ArticleGoogle Scholar
  70. 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.View ArticleGoogle Scholar
  71. Moeser GM, Carrington E. Seasonal variation in mussel byssal thread mechanics. J Exp Biol. 2006;209(10):1996–2003.View ArticlePubMedGoogle Scholar
  72. 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.View ArticleGoogle Scholar
  73. Merkel JR, Dreisbach JH, Ziegler HB. Collagenolytic activity of some marine bacteria. Appl Microbiol. 1975;29(2):145–51.PubMed CentralPubMedGoogle Scholar
  74. 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.Google Scholar
  75. Lima FP, Wethey DS. Three decades of high-resolution coastal sea surface temperatures reveal more than warming. Nat Commun. 2012;3:704.View ArticlePubMedGoogle Scholar
  76. 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.Google Scholar
  77. 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.View ArticleGoogle Scholar


© Zardi et al.; licensee BioMed Central. 2015
