Front. Mar. Sci. Frontiers in Marine Science Front. Mar. Sci. 2296-7745 Frontiers Media S.A. 10.3389/fmars.2017.00075 Marine Science Original Research Congruent Molecular and Morphological Diversity of Macaronesian Limpets: Insights into eco-evolutionary Forces and Tools for Conservation Carreira Gilberto P. 1 * Shaw Paul W. 2 Gonçalves João M. 3 4 McKeown Niall J. 2 1Direção de Serviços de Biodiversidade e Política do Mar, Direção Regional dos Assuntos do Mar, Secretaria Regional do Mar, Ciência e Tecnologia Colónia Alemã, Portugal 2Institute of Biological, Environmental and Rural Sciences, Aberystwyth University Aberystwyth, UK 3MARE – Marine and Environmental Sciences Centre - Azores Horta, Portugal 4OKEANOS-Research Center, University of the Azores Horta, Portugal

Edited by: Tony J. Pitcher, University of British Columbia, Canada

Reviewed by: Sarah Helyar, Queen's University Belfast, Ireland; Ilaria Coscia, University of Salford, UK

*Correspondence: Gilberto P. Carreira gilberto.mp.carreira@azores.gov.pt

This article was submitted to Deep-Sea Environments and Ecology, a section of the journal Frontiers in Marine Science

 16 03 2017 2017 4 75 30 12 2016 02 03 2017 Copyright © 2017 Carreira, Shaw, Gonçalves and McKeown. 2017 Carreira, Shaw, Gonçalves and McKeown

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Combined analysis of genetic and morphological variation can provide synergistic insights into eco-evolutionary forces shaping biodiversity, as well as tools for conservation and management. Macaronesian limpets are undergoing severe declines due to overexploitation which calls for an evaluation of the evolutionary significance and taxonomic status of populations. This study details the analysis of genetic (mtDNA sequencing) and morphological (geometric analysis of shell shape) variation among Macaronesian populations of Patella aspera, P. rustica, and P. candei. In the case of P. aspera and P. rustica, this also included analysis of mainland conspecifics. mtDNA revealed distinct phylogeographic patterns for the three species. For both P. aspera and P. rustica Macaronesian and continental samples were reciprocally monophyletic, with shallower, and distinct, phylogeographic structure within both clades. For P. candei, the Macaronesian endemic, three reciprocally monophyletic groups corresponding to (i) the Azores, (ii) Madeira and (iii) Selvagens and Canary Islands were resolved. The different patterns for each species are compatible with independent processes of colonization and demographic processes. Morphological differentiation matched the major phylogeographic groupings (i.e., Macaronesian v Continent for P. aspera and P. rustica; Azores/Madeira/Selvagens and Canary Islands for P. candei) compatible with such variation being shaped by genetic drift associated with distinct demographic histories of lineages as well as historical/recurrent selection pressures. The genetic/morphological congruence supports the formal recognition of at least three described P. candei subspecies. Moreover, the data prompt recommendation of the separate management of limpets from each archipelago. At present management efforts are hampered by illegal harvesting of limpets. From an applied viewpoint, this study confirms that morphology retains useful information on genetic status and thus represents a potentially cheap method applicable to limpet conservation.

 Gastropods phylogeography evolution conservation phenotype biodiversity Introduction

Patterns of variation in morphological and genetic diversity for a species reflect to varying degrees, historical and recurrent evolutionary processes, and can provide insights into the interplay between neutral and adaptive forces. Such information is fundamental to understand eco-evolutionary processes (Slatkin, 1987) and can be useful for the conservation of intraspecific biodiversity in the case of vulnerable or highly exploited species (Hilborn et al., 2003; Mariani et al., 2012). While advances in molecular biology are providing increased resolution of biocomplexity among marine taxa, there remains a major demand to obtain reliable and informative biological information for a plethora of populations/species at low operational costs for management and conservation purposes (Frankham, 2010). As such, comparative studies of morphological and genetic divergence can still have an important contribution to biodiversity conservation (Strange et al., 2008; Mariani et al., 2012).

The biogeographic region of Macaronesia is comprised of the four ocean archipelagos of the Azores, Madeira, Selvagens and Canaries (the archipelago of Cabo Verde is also sometimes included into this definition). These islands are all the result of volcanic activity and have never been connected to continents, thus their biota are believed to have resulted from colonization and in situ evolution. Across that region, Patellid limpets represent an important socioeconomic resource and have likely been exploited since the islands were first colonized (Hawkins et al., 2000), which, in the case of Canaries, goes back to around the middle of the 1st millennium B.C.E with recognized effects on the abundance of P. candei, at least in one island of the archipelago (Gonzalez-Lorenzo et al., 2015) (the human colonization of the Azores and Madeira goes back to the 15th century). The high levels of exploitation have contributed to severe declines in limpet abundance with many populations considered to be overexploited (Martins et al., 2008). The reductions in limpet abundance is also reported to be altering ecosystem services of the coastal areas, as it has been linked to the increase in turf-forming algae (Martins et al., 2010).

Three species of Patella are recognized on the Macaronesian islands: Patella aspera, P. rustica, and P. candei (Christiaens, 1973; Hawkins et al., 1990; Ridgway et al., 1998; Koufopanou et al., 1999). P. aspera is the widest distributed Patella species. It inhabits the lower intertidal and sub-tidal down to 6 m deep and is distributed across most of Europe, Mediterranean, North Africa, as well as the Azores, Madeira, Selvagens, and Canaries. Continental populations have been referred to as P. ulissyponensis, while Macaronesian populations have been referred to as P. aspera (Hawkins et al., 2000; Weber and Hawkins, 2005), sometimes also called P. u. aspera (Ferraz et al., 2001). P. rustica occurs throughout Europe, up to northern Iberia, including the Mediterranean and North West African Atlantic coasts south to Mauritania. P. rustica inhabiting the Macaronesian Islands are sometimes referred to as P. piperata (Christiaens, 1968; Côrte-Real et al., 1996b). This species is absent from the Azores and is believed to have undergone a recent range expansion across the western coast of Iberia due to climate change (Lima et al., 2006). P. candei is endemic to the Macaronesian Islands, where it inhabits the intertidal fringe, usually less than 3 m depth. This species overlaps in its habitat with other species, according to the archipelago. In Madeira and Canaries, two species share partially the same habitat with P. candei: P. rustica (upper intertidal) and P. aspera (low intertidal and subtidal). In the Azores, however, only P. aspera partially shares the same habitat with P. candei, in the low intertidal and upper subtidal, as P. rustica is absent from that archipelago.

Phylogeographic patterns among Macaronesian limpets are well established (Sá-Pinto et al., 2008) and reveal ancestral splits and contemporary restricted gene flow from mainland forms (for P. aspera and P. rustica), as well as further restricted gene flow among archipelagos (for all three species). In contrast, informative quantitative morphological data are lacking, largely due to difficulties in identifying appropriate and informative landmarks on shells (e.g., Cabral and Jorge, 2007). Where morphology reflects genetic patterns, preliminary morphological surveys can motivate and inform future studies into factors shaping gene flow. Furthermore, in such cases phenotypic traits can be used as valid criteria to establish management units (Garnier et al., 2005). Additionally, morphological differentiation can reveal cryptic patterns of incipient genetic divergence (e.g., that may not be detected at assayed genetic markers that have not attained migration-drift equilibrium or sufficient statistical power (Nice and Shapiro, 1999). Generation of robust morphological data is particularly relevant to the taxonomy of P. candei for which early morphological patterns prompted the description of four geographically distinct subspecies (P. c. gomesi in Azores; P. c. ordinaria in Madeira, P. c. candei in Selvagens, and P. c. crenata in Canaries (Christiaens, 1973). However, subsequent studies have disagreed with various aspects of this taxonomy and it remains controversial (Côrte-Real et al., 1996a,b; Weber and Hawkins, 2002; Sá-Pinto et al., 2008).

This study details the analysis of genetic (mtDNA sequencing) and morphological variation among Macaronesian populations of P. aspera, P. rustica and P. candei. In the case of P. aspera and P. rustica, this also included analysis of mainland conspecifics. Morphological variation was assessed using the first, to our knowledge, implementation of a geometric approach for limpets and particular attention was paid to the morphological diagnosability of members of different lineages. This information is then used to discuss (i) the factors shaping biodiversity and (ii) the potential utility of such an inexpensive morphological approach as a management tool. The patterns are also interpreted in the context of P. candei taxonomy.

 Materials and methods

Specimens of P. aspera, P. rustica, and P. candei were collected from across their distribution ranges across the Azores, Madeira, Selvagens and Canaries. For P. rustica and P aspera, representative continental samples were also collected and examined (Table 1). Samples were collected randomly by hand, throughout the intertidal rocky shore, during low tide periods. Tissues were preserved in ethanol (96%) and shells cataloged for further morphometric analysis. In the case of shells, only non-broken and not severely eroded specimens were considered for morphometric study.

Number of specimens examined for the three Patella species, localities of collection and data (COI and cytochrome b, or CYTB) and shell shape (2-D conformations, composed of 27 landmarks).

Species Geographical distribution Area/Archipelago Locality/Island COI CYTB COI-CYTB Morphology
Patella candei Azores Azores Corvo 25 27 15 43
Madeira Flores 8 19 4 26
Selvagens Graciosa 11 8 4 20
Canary Islands São Jorge 13 11 4 58
Terceira 13 16 5 38
Faial 12 16 2 55
Pico 17 19 4 206
São Miguel 10 16 4 89
Santa Maria 10 12 – 57
Formigas Islets 14 18 9 §
Madeira Madeira 14 17 10 126
Porto Santo 21 16 3 38
Deserta 19 17 10 98
Selvagens Selvagem Grande 20 20 6 §
Canary Islands Gran Canária 22 19 13 34
Lanzarote - Punta Pechiguera 13 23 – 35
Lanzarote - El Cochino, Tymanfaya 14 14 10 21
La Graciosa 5 14 3 41
P. rustica Madeira, Madeira Deserta 3 3 3 4
Selvagens, Porto Santo 10 14 4 21
Canary Islands, Canary Islands Gran Canária 14 10 8 18
North west Africa, Lanzarote 10 19 6 16
Atlantic SW Europe, La Graciosa 12 8 6 14
West Mediterranean, North west Africa Mauritania 1 3 3 4
Aegean Sea. Atlantic SW Europe, Algarve 16 16 15 15
West Mediterranean Ceuta 14 23 13 44
Tuscany 5 13 4 18
Aegean Sea Greece 19 13 10 30
P. aspera Azores, Azores Corvo 16 9 8 43
Madeira, Flores 2 2 2 §
Selvagens, São Jorge 26 13 12 §
Canary Islands, Faial 13 13 10 15
Atlantic Europe, Terceira 8 5 2 29
West Mediterranean, Pico 25 18 13 20
Aegean Sea. São Miguel 19 9 5 28
Madeira Madeira 12 4 1 59
Porto Santo 14 11 5 17
Deserta 4 4 4 2
Selvagens Selvagem Grande 11 5 4 4
Canary Islands Gran Canária 13 18 4 37
La Graciosa 1 3 – 2
Lanzarote 14 2 1 15
Atlantic Europe Algarve 10 21 6 13
Costa da Caparica 12 17 9 42
Figueira da Foz 15 12 12 13
Devon, Great Britain 1* – – –
West Mediterranean Ceuta 5 3 3 §
Tuscany 7 10 5 4
Aegean Sea Greece 11 13 5 19
Total 611 606 289 1564

Tissue sample obtained from a voucher deposited in the Natural History Museum of London, UK; § locations for which the shells were not included for analysis, due to the fact that these were highly eroded shells.

 Genetic analysis

Where possible sequences were obtained for each individual for the cytochrome oxidase I gene (~650 bp fragment PCR amplified using universal primers LCO1490 and HCO2198, Folmer et al., 1994) and cytochrome B gene (~350 bp PCR amplified using degenerate primers UCYTB151F and UCYTB270R, Merritt et al., 1998). PCR products were cleaned using ExoSap and sequenced in both directions on an Applied Biosystems sequencer using Big Dye technology. Sequences were trimmed and edited manually in CHROMAS (http://www.technelysium.com.au) and the alignments were analyzed with the programs ESEE (Cabot and Beckenbach, 1989), BioEdit 7.0.9.0 (Hall, 1999), DnaSP 4.50.3 (Rozas et al., 2003), MEGA 4.0.2 (Tamura et al., 2007), and DAMBE 5.0.52 (copyright Xuhua Xia, http://dambe.bio.uottawa.ca/Include/software.aspx). Neighbor joining P-uncorrected phylogenies were inferred by using PAUP* Neighbor joining P-uncorrected trees obtained in MEGA. The selection of the best DNA substitution model was performed for the concatenated datasets through MrMODELTEST 2.3 (Nylander, 2004), jointly with PAUP* 10b (Swofford, 2003). Maximum likelihood (ML) phylogenies for the COI-CYTB composite (concatenated sequence between COI and cytochrome b) haplotypes were inferred with the routine implemented in PhyML, as described by Guindon and Gascuel (2003). ML phylogenies were assessed by bootstrapping, with 2,000 pseudo replicates (Felsenstein, 1985). The ML trees were mid-point rooted, visualized and edited using FigTree 1.2.2 (Andrew Rambaut, 2006–2009, http://tree.bio.ed.ac.uk/). Genetic variation was described using indices of haplotype and nucleotide diversity (Nei and Tajima, 1981; Nei, 1987; h and ∏, respectively), extracted through DnaSP.

 Morphological analysis

Each shell was photographed with a digital camera (Nikon CoolPix 775, with lens Nikkor 5.8–7.4 mm, I:2.8–4.9) fixed to a tripod. Each shell was superimposed to a fan composed of a system of axes equally distant in 15°, and centered to a frame of size 110 × 70 mm (Figure 1). Each shell was orientated with the anterior end to the left and two additional lines were used on the right side, distant to the horizontal line in ±22°. TPS data files were prepared with tpsUtils (Rohlf, 2006a) and the digitalization of each specimen was made with tpsDig (Rohlf, 2006b). This system allowed the collection of 27 landmarks in each shell (one located on the shell apex, and the remaining landmarks located throughout the shell border). The pictures were individually scaled and a dataset composed of shape conformations obtained, consisting of two-dimensional (x,y) coordinates of the 27 landmarks that were collected from each specimen (Figure 1).

Scheme for shape data collection on limpet shells considering 27 landmarks. In this case, a specimen of P. candei collected on Pico Island (Azores).

The individual conformations were aligned with the generalized Procrustes superimposition method, as implemented in tpsRelW (Rohlf, 2006c). In order to relax the assumption of homology among landmarks, the border points were allowed to slide as semi-landmarks along the shell border, according to a pre-established sliders file, prepared for each dataset with tpsUtils. A matrix of shape descriptors generated on tpsRelw (Rohlf, 2006c), known as the weight matrix, was extracted with tpsRelw. This matrix is composed of the partial warps (PW) and the two uniform (affine) shape components. First, a relative warp analysis (RWA) was performed. For this analysis, and to avoid the noise provided by the typical (and possibly random) reticulated nature of the shell border in limpets, we set α parameter of the RWA to 1, which drives the analysis to focus primarily on the overall information on shape, when performing a RWA. The matrix of relative warps (RW) was analyzed by plotting means and 95% confidence intervals for each group previously defined as a molecular lineage.

One-way analyses of variance (ANOVA) between groups (mitochondrial lineages) were performed and Bonferroni pairwise tests were computed in order to statistically assess shape differences between pairwise groups, as described by the first three RW components. A visual representation of the shape variation described by RW was obtained by plotting thin-plate-spline deformation grids associated with extreme values of RW for each group, which consist of a graphical illustration of the bending energy matrix, as described by each RW component.

 Descriminant analysis and assessment of molecular lineage/morphological congruence

The weight matrixes were imported to the STATISTICA®—data analysis software system, version 6.0 (http://statistica.software.informer.com/6.0/, Inc. 2001), and a generalized discriminant analysis performed by defining the groups according to the mitochondrial lineages. In addition, cross-validation procedures were performed by previously keeping 25% (P. candei) or 33% (P. rustica and P. aspera) (taken randomly from the specimen groups included in each clade) as a training sample.

In order to visually inspect the shape variability associated with the inferred discriminant model, multivariate regression analyses were computed of the weight matrix (set as dependent variable matrix) on the coefficient of variate matrix (CV, set as the matrix of independent variables). The fit of these regression analyses were assessed through the significance of the multivariate Goodall F, based on a re-sampling procedure (5,000 permutations). These analyses were performed on tpsRegr (Rohlf, 2006d), which also allowed the visualization of the correspondent deformation grids, under a thin-plate-spline model, based on the bending energy matrix. The congruence between molecular lineage/morphology was assessed by evaluating the degree of molecular disjunction between clades and the statistical grouping both for RW and for discriminant analysis (cross-validation).

 Results Patella aspera

In 240 individuals of P. aspera that were screened for COI, 46 haplotypes were detected (h = 0.82 ± 0.020), and 34 haplotypes were detected in 175 individuals that were screened for cytochrome b (h = 0.82 ± 0.020). The concatenated dataset was composed of 39 different haplotypes. The nucleotide diversity of COI for the complete alignment was 0.02165 ± 0.00124. The nucleotide diversity of the complete cytochrome b alignment was 0.03350 ± 0.00128.

Both the neighbor joining p-uncorrected phylograms and the ML inferred phylogeny revealed two reciprocally monophyletic groups, one composed of the haplotypes sampled across the Macaronesian islands and the second containing haplotypes sampled across the Atlantic and Mediterranean continental coasts (Figure 2). Within both clades there was further phylogeographic structuring. Among the Macaronesian samples, the Azores shared no haplotypes with the other archipelagoes, while there was also no haplotype sharing between the Atlantic and Mediterranean samples. Canaries and Madeira shared haplotypes for both COI and cytochrome b.

Neighbor joining P-uncorrected phylograms for COI and CYTB haplotypes, and ML tree among concatenated haplotypes in P. aspera. Values located next to the branches indicate percentage of bootstrap replicates that supported each branch (2,000 bootstrap replicates).

Patterns of shell shape variation showed concordance with the inferred phylogeographic patterns: the main divergence in terms of shell shape was found between the two main mitochondrial phylogenetic lineages (Macaronesia/continent), particularly demonstrated for the RW1 and the uniform (affine) shape descriptors (Figure 3). The ANOVA clearly split the insular from the Atlantic and Mediterranean continental populations, in terms of polygonality of the shells (ANOVA, RW1, F = 48.26248, df = 358, P = 0.000013) as well as the lateral vs. longitudinal compression (ANOVA, Ux, F = 48.46422, df = 358, P = 0.000012). The RW1 characterized the polygonality, typical of the Macaronesian shells, by the eccentricity of the landmarks 9, 13, 17, 26, and 27, as well as a more posteriorly positioned apex. The continental shells, on the other hand, presented a smoother and less angular shell shape. The affine shape descriptors showed the continental shells as laterally compressed, in opposition to the Macaronesian shells, which were wider and longitudinally compressed. This differentiation was evident in the high percentage of correct assignment of individuals to either the Macaronesian (83.3%) or continental (80.0%) groups (overall 84.4% correct assignment to regional group) (Figure 4, Table 2). Lower shell shape differentiation was reported among the lineages within the Macaronesian and continental groups. A gradient of shape was obtained (Figure 4) ranging from the continental shells, characterized by only slight polygonal and laterally compressed shell shapes, to the Macaronesian shells, characterized by highly polygonal and laterally wider/longitudinally compressed shell shapes.

Mean and 95% confidence interval of the relative warp scores and first uniform (affine) shape component, for individuals from the four identified P. aspera mitochondrial lineages. The linking lines indicate non-significant differences, according to the pairwise Bonferroni post-hoc multiple comparison (p > 0.05). Scatter plots of the canonical variates (roots 1, 2, and 3) obtained by discriminant analysis of the shape descriptors (weight matrix) of P. aspera, and considering the four mitochondrial lineage groups. Goodall F = 30.1909, df = 150, 17,900, P = 0.0000.

Canonical variates obtained from discriminant analysis, when analyzing the two main clades of P. aspera: Atlantic Islands (composed of the two reciprocally monophyletic groups of Madeira and Canary Islands); and Continent (comprising the sample localities sampled across the Atlantic European, Atlantic African, Mediterranean African and Mediterranean European coasts). A cline on the shell apex position is visible, from more anterior (Atlantic Islands) to more central (Continent). Goodall F = 62.7807, df = 50, 18,000, P = 0.0000.

Cross-validation reassignments, based on the GDA of the weight matrixes, composed of the non-uniform (PW) and uniform (affine) shape descriptors for each species.

  P. aspera - Correct cross-validation (train sample = 33%), Canonical R = 0.77343, Eigen-value = 1.487927, Wilk's lambda = 0.186704, χ2 = 359.1412, df = 150, P = 0.00000.
Lineage Mediterranean Atlantic Azores Madeira and Canary Islands Total (%) Main disjunction (%)
Mediterranean 2 3 1 1 29 84
Atlantic 9 10 2 2 43
Azores 1 6 25 13 56 80
Madeira and Canary Islands 3 4 13 25 56
Total 15 23 41 41 52
  P. rustica - Correct cross-validation (train sample = 33%), Canonical R = 0.869988, Eigen-value = 3.113166, Wilk's lambda = 0.068103, χ2 = 255.2398, df = 150, P = 0.00000.
Lineage Atlantic and west Mediterranean Eastern Mediterranean Madeira Canary Islands Total (%) Main disjunction (%)
Atlantic and west Mediterranean 16 5 3 3 59 84
Eastern Mediterranean 5 5 0 0 50
Madeira 0 1 4 3 50 83
Canary Islands 0 3 6 7 43
Total 21 14 13 13
  P. candei - Correct cross-validation (train sample = 25%), Canonical R = 0.8603; Eigen-value = 2.8483; Wilk's lambda = 0.1481; χ2 = 1286.3080; dF = 100; P = 0.00000.
Lineage Azores Madeira Canary Islands – Total(%) –
Azores 132 0 3 – 98 –
Madeira 4 55 6 – 85 –
Canary Islands 5 3 25 – 76 –
Total 141 58 34 – 91 –
P. rustica

In 110 specimens of P. rustica that were screened for COI, 41 haplotypes were found (h = 0.94 ± 0.011), and 30 cytochrome b haplotypes were found in 115 specimens (h = 0.83 ± 0.024). 39 contiguous haplotypes were identified. Inferred phylogenies revealed a hierarchically clustering into four monophyletic phylogeographic groups (Figure 5). The major phylogenetic split separated the Macaronesian and continental lineages. Within both clades reciprocally monophyletic units composed of (i) the Madeira and Canary islands and (ii) Aegean Sea and Atlantic/west Mediterranean were also resolved.

Neighbor joining P-uncorrected phylograms for COI and CYTB haplotypes, and ML tree among concatenated haplotypes in P. rustica. Values located next to the branches indicate percentage of bootstrap replicates that supported each branch (2,000 bootstrap replicates).

Shell shape variation for P. rustica also aligned with the phylogeographic patterns in clearly separating the Macaronesian and continental samples, whilst also reporting a lower level of diagnosability of members of the sub-clades (i.e., Madeira, Canary Islands, Atlantic/west Mediterranean, Aegean Sea). The RW1 (which explained 48.15% of the variation in shell shape) separated Macacronesian and continental shells as well as shells from the Canary Islands and Madeira, (Figure 6). The main character that was responsible for this observed pattern was the relative location of the shell apex, which was more anterior among Madeira individuals, and to a lesser degree in Canary Island individuals. The shells from the Atlantic and Mediterranean showed a more relatively centered apex location. This pattern was also apparent in RW3 (19.10% of shell shape variability). The shells from Madeira were also significantly more laterally compressed, when compared to the other three lineages, which did not significantly differ (P > 0.05) by presenting average values of the uniform (affine) shape component.

Above: Mean and 95% confidence interval of the relative warp scores and first uniform (affine) shape component, for each P. rustica mitochondrial lineage. The linking lines indicate non-significant differences, according to the pairwise Bonferroni pos-hoc multiple comparison (p > 0.05). Deformation grids were based on the matrix of bending energy, under the model of thin-plate spline, and correspond to the most extreme axe RW scores. Below: Scatter plots of the canonical variates (roots 1, 2, and 3) obtained by discriminant analysis. Deformation grids obtained by regression (F = 52.4915, df = 150, 9,000, P = 0.0000).

The four lineages were also separated significantly (Canonical R = 0.862436, Eigen-value = 2.903149, Wilk's lambda = 0.109511, χ2 = 345.0298, df = 150, P < 0.00001) based on the discriminant model inferred by GDA, although a pattern of partial overlap was found (Figure 6). The first canonical root (73%) clearly separated Macaronesian and continental shells, based on the more anterior position of the apex. The second canonical root, on the other hand, separated shells from the two island populations, with the ones from the Madeira archipelago presenting a more anterior apex and being slightly laterally compressed. The third canonical root (9%), although showing some overlap, tended to separate shells based on the anterior border, which was anteriorly pointed in Aegean Sea shells compared to the Atlantic and Western Mediterranean. Polygonality was not apparent on P. rustica shells.

The shape overlap found in both the RWA and GDA was reflected in the results of the cross-validation assignment tests wherein high levels of assignment success were obtained between the Macaronesian and continental samples but assignment success was lower when classification was to the four mtDNA clades i.e., Madeira/ Canary islands/Atlantic/Mediterranean) (Table 2). The GDA, when performed by assuming these two main mitochondrial lineages as discriminant groups, was significant (Goodall F = 127.8290, df = 50, 9100, P < 0.0001) and showed a clear separation (Figure 7) based on the canonical root. The thin-plate-spline regression of the weight matrix on the canonical root confirmed the apex position as the main character separating these two lineages: Macaronesian specimens presented significantly anterior apex locations when compared to the continental populations.

Plot of the canonical variates obtained with discriminant analysis, when analyzing the two main clades of P. rustica: Atlantic Islands (composed of the two reciprocally monophyletic groups of Madeira and Canary Islands); and Continent (comprising the sample localities sampled across the Atlantic European, Atlantic African, Mediterranean African, and Mediterranean European coasts). A cline on the shell apex position is visible, from more anterior (Atlantic Islands) to more centrally located (Continent). Goodall F = 127.8290, df = 50, 9100, P = 0.0000.

 P. candei

In 261 individuals of P. candei that were screened for COI (628 bp), 81 haplotypes were detected (h = 0.93 ± 0.010). For cytochrome b (308 bp) 40 haplotypes were detected in 306 screened individuals of P. candei (h = 0.73 ± 0.022).The nucleotide diversity (π) in the Azores was 0.00548 ± 0.00051 for this species and for cytochrome b was π = 0.04822 ± 0.00373. The concatenated dataset was composed of 58 composite haplotypes. Phylogenetic analysis revealed a strong partitioning of three phylogeographic groups corresponding to (i) the Azores, (ii) Madeira and (iii) Selvagens and Canary Islands (Figure 8).

Neighbor joining P-uncorrected phylograms for COI and CYTB haplotypes, and ML tree among concatenated haplotypes in P. candei. Values located next to the branches indicate percentage of bootstrap replicates that supported each branch (2,000 bootstrap replicates).

The complete disjunction that was detected among the three geographically-based mitochondrial lineages of P. candei was also resolved in the shell shape (Figure 9). The relative warp analysis completely separated the three lineages, based on the polygonal degree (RW1, F = 24.783, DF = 931, P = 0.0000), on the eccentricity of the apex position (RW3, F = 29.420, DF = 931, P = 0.0000), and on the shell lateral compression (Ux, F = 21.818, DF = 931, P = 0.0000). Only RW2 displayed non-significant differences among the three mitochondrial lineages (ANOVA, F = 2.874, DF = 931, P = 0.0570) (Figure 9).

Relative warp scores and first uniform (affine) shape component (left), for each P. candei mitochondrial lineage. The linking lines indicate non-significant differences, according to the pairwise Bonferroni pos-hoc multiple comparison (P > 0.05). Deformation grids were based on the matrix of bending energy, under the model of thin-plate spline, and correspond to the most extreme axe RW scores. The grids for the affine component were obtained by setting α = 0, when performing the relative warp analysis. Scatter plot of the canonical variates (roots) obtained by discriminant analysis of the shape descriptors (weight matrix) of P. candei (right), and considering the mitochondrial lineage groups. Deformation grids obtained through multivariate regression. Goodall F = 67.2600, DF = 100, 46550, P < 0.0001.

This morphometric analysis revealed that the lineages inhabiting the Madeira and Canary archipelagos are typified by clear polygonal shell shapes, and that a relatively anterior positioned apex characterizes individuals in the lineage inhabiting the Azores. The affine shape descriptor showed a gradient of lateral compression, from Canary Islands (more rounded and longitudinally compressed) to the Azores where the shells tended to be more laterally compressed.

The discriminant analysis performed among molecular lineages was significant (Canonical R = 0.8572, Eigen-value = 2.7714, Wilk's lambda = 0.1520, χ2 = 1707.667, DF = 100, P < 0.0001), and the canonical variate roots separated the three lineages based on the following shape characters (Figure 9): the first canonical root (which described 78.8% of the shape variation) primarily separated the Azores from Madeira, based on the degree of polygonality. Shells from Madeira presented highly polygonal shells, characterized by the eccentricity of the landmarks 2, 9, 17, 23, 26, and 27. A number of landmarks showed concentric displacements, such as landmarks 1, 3, and 8. The shells from the Azores lineage, on the other hand, showed a more oval general shape. The shells from the Canary Islands were positioned intermediate between the other two lineages, and were positively separated based on the second canonical root (which accounted for 21.2% of the shape variation), characterized by a posteriorly positioned apex and a more reticulated shell border. Landmark 13 was typically positioned eccentrically.

The cross-validation percentage of correct shell-shape reassignment of specimens to their true phylogenetic lineage was high (Table 2), especially for Azores individuals for which only 3 specimens (2% of the training sample) were misclassified as belonging to the Canary Islands clade.

 Discussion

Macaronesian limpets are undergoing severe declines due to overexploitation which calls for an evaluation of the evolutionary significance and taxonomic status of populations. This study reports considerable evolutionary divergence of the Macaronesian limpets, revealed through highly concordant patterns of genetic and morphological variation. mtDNA revealed distinct phylogeographic patterns for the three species compatible with historical and contemporary gene flow restrictions previously reported by Sá-Pinto et al. (2008). Specifically, for both P. aspera and P. rustica Macaronesian and continental samples were reciprocally monophyletic, with lower levels of phylogeographic structure occurring within both clades. For P. candei, the Macaronesian endemic, three reciprocally monophyletic groups corresponding to (i) the Azores, (ii) Madeira and (iii) Selvagens, and Canary Islands were resolved. The different patterns found among these species are compatible with independent processes of colonization and demographic processes for the three species.

Geometric morphometric analysis detected significant morphological differentiation that matched major phylogeographic groupings (i.e., Macaronesian v Continent for P. aspera and P. rustica; Azores/Madeira/Selvagens and Canary Islands for P. candei). While plasticity is described for limpets (Hawkins et al., 1990) and widely reported for molluscs (Trussell, 2000) such concordant genetic and morphological differentiation is often taken as indicating the roles of genetic drift and/or natural selection. While a greater understanding of the relative roles of selection and plasticity could be obtained through comparison of QST and FST values (but see Whitlock, 2008) this requires a common garden approach. We suggest that the results strongly indicate that morphological variation has in some way been shaped by genetic drift associated with distinct demographic histories of lineages as well as historical/recurrent selection pressures.

The development of geometric morphometric techniques has considerably advanced morphological studies (Davis et al., 2016) and in comparison to traditional distance-based methods permits a more detailed description of phenotypic variation. For both P. aspera and P. rustica, Macaronesian and continental populations were observed to be highly genetically diverged and reciprocally monophyletic. These major genetic breaks coincided with the largest components of morphological divergence. For P. aspera the continental shells were characterized as being less polygonal and more laterally compressed than their Macaronesian counterparts which showed a trend to become more polygonal and reticulated, presenting a lateral expansion. For P. rustica the relative apex position appeared to be the main character difference between continental and Macaronesian shells. Among the Archipelagos both species reported differing levels of genetic distinctiveness. For P. aspera there were clear haplotype frequency differences among sites, but not reciprocal monophyly, whereas for P. rustica there was a congruent complete lineage sorting between Madeira and Canary islands. These differing depths of divergence were reflected in the observed morphological differentiation. For P. aspera individuals could not be morphologically discriminated across islands, whereas for P. rustica morphological differentiation between Madeira and Canaries was significant, though lower than that reported in comparisons with continental conspecifics, which also reflects the hierarchical patterns resolved by mtDNA.

Despite the general correspondence between genetic and morphological data, P. rustica revealed an interesting incongruence. Whereas morphological differentiation was detected between the genetically distinct archipelagos no morphological differentiation was detected between Aegean Sea and the Western Mediterranean /Atlantic lineages despite the comparably high genetic divergence. While it must be stressed that only one sampling location from the Aegean Sea has been included in this study, such a morphological similarity could reflect stasis due to stabilizing selection. Additionally, this result could reflect the accelerated evolution of Macaronesian limpet diversity due to historical demographics and selection factors.

P. candei exhibited the highest levels of clade divergence, with three reciprocally monophyletic groups. Sá-Pinto et al. (2008) interpreted such divergence as indicating single colonization events, and subsequent genetic isolation for each archipelago. Such isolation, along with post colonization adaptation, would be expected to drive phenotypic divergence and this pattern was supported by morphometrics which reported the three lineages to exhibit a greater level of morphological differentiation, than observed for the other species. The shells carried by limpets inhabiting Madeira are characteristically polygonal, when compared with the shells of the limpets that inhabit the Azores. The shells from the Canary Islands, on the other hand, present a typical reticulated border and with its shell apex centrally located, when compared with the other two lineages of P. candei. These results confirm the high morphological differentiation described (see for example Christiaens, 1973) between the three subspecies P. c. gomesi (Azores), P. c. ordinaria (Madeira) and P. c. crenata (Canary Islands) of P. candei. Specimens from P. c. candei (Selvagens) were not included in the analysis here, given their eroded status, however, subsequent studies suggest this eroded phenotype may have a genetic component (Carreira, 2010). While mtDNA has failed to differentiate P. candei from Canaries and Selvagens, this eroded phenotype may reflect reproductive isolation and incipient divergence that may be detectable with more powerful genetic markers such as microsatellites.

Islands are regarded as excellent model systems to investigate the relative roles of historical contingency and ecological determinism in shaping biodiversity (Emerson, 1985). Cunha et al. (2008) reported replicated biogeographic patterns between large and small coned venomous snails in Cabo Verde wherein ecological determinism, specifically recurrent gene flow restrictions and adaptation, had overcome the effects of different evolutionary histories. In this study, the most pronounced phylogeographic and morphological diversity among archipelagos was reported for P. candei, with lower levels of genetic/morphological divergence reported for the other species. While the distinct phylogeographic structures support unique colonization histories for the three species (Sá-Pinto et al., 2008), as P. candei are the oldest of the three species on Macaronesia, it may be that the differing patterns among the species in Macaronesia reflect different temporal stages of divergence by post-colonization determinism toward a replicated, albeit general, model of among archipelago diversification. As such the Macaronesian limpets represent an excellent model system for studying the stages of evolutionary divergence.

A fundamental aim of conservation is to identify and prioritize evolutionary significant components of biodiversity. Comparative genetic and morphometric analysis here highlight the unique genetic and morphological characteristics of the Macaronesian limpets, and implicate underpinning roles for historical contingency and contemporary determinism in the evolution of this biodiversity. The subtlety of evolutionary processes can complicate conservation efforts focused at the sub-species rank (Crandall et al., 2000), however the genetic and morphological divergence support the formal recognition of at least three described P. candei subspecies. For all three species the genetic divergence among archipelagos must be considered a conservative reflection of contemporary isolation, and it is recommended that each archipelago be managed separately as such isolation increases the vulnerability of populations, which in the case of P. candei represents a real extinction risk. At present management efforts are severely hampered by illegal harvesting of limpets. From an applied viewpoint, this study confirms that morphology retains useful information on genetic status and thus represents a potentially cheap method applicable to limpet forensics.

 Data accessibility

All sequences have been deposited on GenBank (accession numbers: KY674542-KY674813).

 Author contributions

All the four authors met the four criteria of authorship, namely: Substantial contributions to the conception or design of the work (GC, PS, and JG); or the acquisition (GC), analysis (GC and NM), or interpretation of data for the work (GC and NM); Drafting the work or revising it critically for important intellectual content (All); Final approval of the version to be published (All); Agreement to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved (All).

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This study was conducted at Royal Holloway, University of London and at the Departamento de Oceanografia e Pescas, da Universidade dos Açores. This study benefited from the scholarship SFRH/BD/6897/2001, through the Portuguese Foundation for Science and Technology and the following grants: OGAMP (Interreg IIIB-MAC/4.2/A2), MARMAC (Interreg 03/MAC/4.2/A1) and BANGEN PCT/MAC/1/C070—Banco genético marinho da Macaronésia. All the samples were collected according to the law enforced on each locality at the time of collection. In Canary Islands, samples were collected during the survey campaigns of the project “Canarias por una costa viva,” in collaboration with the team of the University of Las Palmas (Gran Canaria). Finally, the open access publication was supported by Portuguese “Fundação para a Ciência e a Tecnologia (FCT)” through the strategic project (FCT/UID/MAR/04292/2013) granted to MARE.

 References Cabot E. L. Beckenbach A. T. (1989). Simultaneous editing of multiple nucleic-acid and protein sequences with ESEE. Comp. Appl. Biosci. 5, 233–234. 10.1093/bioinformatics/5.3.2332766009 Cabral J. P. Jorge R. M. N. (2007). Compressibility and shell failure in the European Atlantic Patella limpets. Mar. Biol. 150, 585–597. 10.1007/s00227-006-0379-0 Carreira G. P. (2010). Molecular Genetic Diversity and Shell Shape Variation in Patella limpets (Mollusca: Patellogastropoda): Evolutionary Inferences and Tools for Species Conservation in the North East Atlantic Ocean Archipelagos. PhD. thesis, University of the Azores. Horta. Christiaens J. (1968). Validité du nom Patella piperata Gould. Bulletin du Muséum National D'Histoire Naturelle II-40, 366–373. Christiaens J. (1973). Révision du genre Patella (Mollusca, Gastropoda). Bull. Muséum Natl. D'Histoire Nat. 3, 1305–1391. Côrte-Real H. Hawkins S. J. Thorpe J. P. (1996a). An interpretation of the taxonomic relationship between the limpets Patella rustica and P-piperata. J. Mar. Biol. Assoc. U. K. 76, 717–732. 10.1017/S0025315400031416 Côrte-Real H. Hawkins S. J. Thorpe J. P. (1996b). Population differentiation and taxonomic status of the exploited limpet Patella candei in the Macaronesian islands (Azores, Madeira, Canaries). Mar. Biol. 125, 141–152. 10.1007/BF00350768 Crandall K. A. Bininda-Emonds O. R. Mace G. M. Wayne R. K. (2000). Considering evolutionary processes in conservation biology. Trends Ecol. Evol. 15, 290–295. 10.1016/S0169-5347(00)01876-010856956 Cunha R. L. Tenorio M. J. Afonso C. Castilho R. Zardoya R. (2008). Replaying the tape: recurring biogeographical patterns in Cape Verde Conus after 12 million years. Mol. Ecol. 17, 885–901. 10.1111/j.1365-294X.2007.03618.x18179424 Davis M. A. Douglas M. R. Collyer M. L. Douglas M. E. (2016). Deconstructing a species-complex: geometric morphometric and molecular analyses define species in the western rattlesnake (Crotalus viridis) PLoS ONE 11:e0146166. 10.1371/journal.pone.014616626816132 Emerson B. C. (1985). Evolution on oceanic islands: molecular phylogenetic approaches to understanding pattern and process. Mol. Ecol. 11, 951–96612030975 Felsenstein J. (1985). Confidence-limits on phylogenies - an approach using the bootstrap. Evolution 39, 783–791. 10.2307/2408678 Ferraz R. R. Menezes G. M. Santos R. S. (2001). Limpet (Patella spp.) (Mollusca: Gastropoda) exploitation in the Azores, during the period 1993–1998. Arquipélago-Life Mar. Sci. Part B. (Suppl 2), 59–65. Folmer O. Black M. Hoeh W. Lutz R. Vrijenhoek R. (1994). DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3, 294–299. 7881515 Frankham R. (2010). Challenges and opportunities of genetic approaches to biological conservation. Biol. Conserv. 143, 1919–1927. 10.1016/j.biocon.2010.05.011 Garnier S. Magniez-Jannin F. Rasplus J. Y. Alibert P. (2005). When morphometry meets genetics: inferring the phylogeography of Carabus solieri using Fourier analyses of pronotum and male genitalia. J. Evol. Biol. 18, 269–280. 10.1111/j.1420-9101.2004.00854.x15715833 Gonzalez-Lorenzo G. Hernandez E. M. Perez-Dionis G. Hernandez A. B. Santos B. G. Diez J. B. (2015). Ineffective conservation threatens Patella candei, an endangered limpet endemic to the Macaronesian islands. Biol. Conserv. 192, 428–435. 10.1016/j.biocon.2015.11.005 Guindon S. Gascuel O. (2003). A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol. 52, 696–704. 10.1080/1063515039023552014530136 Hall T. A. (1999). BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids. Symp. Ser. 41, 95–98. Hawkins S. J. Côrte-Real H. Pannacciulli F. G. Weber L. C. Bishop J. D. D. (2000). Thoughts on the ecology and evolution of the intertidal biota of the Azores and other Atlantic islands. Hydrobiologia 440, 3–17. 10.1023/A:1004118220083 Hawkins S. J. H. B. Côrte-Real S. M. Martins H. R. Santos R. S. Frias Martins A. M. (1990). A note on the identity of Patella in the Azores. Açoreana Supplement, 167–173. Hilborn R. Quinn T. P. Schindler D. E. Rogers D. E. (2003). Biocomplexity and fisheries sustainability. Proc. Natl. Acad. Sci. U.S.A. 100, 6564–6568. 10.1073/pnas.103727410012743372 Koufopanou V. Reid D. G. Ridgway S. A. Thomas R. H. (1999). A molecular phylogeny of the patellid limpets (Gastropoda: Patellidae) and its implications for the origins of their antitropical distribution. Mol. Phylogenet. Evol. 11, 138–156. 10.1006/mpev.1998.055710082617 Lima F. P. Queiroz N. Ribeiro P. A. Hawkins S. J. Santos A. M. (2006). Recent changes in the distribution of a marine gastropod, Patella rustica Linnaeus, 1758, and their relationship to unusual climatic events. J. Biogeogr. 33, 812–822. 10.1111/j.1365-2699.2006.01457.x Mariani S. Peijnenburg K. Weetman D. (2012). Independence of neutral and adaptive divergence in a low dispersal marine mollusc. Mar. Ecol. Prog. Ser. 446, 173–187. 10.3354/meps09507 Martins G. M. Jenkins S. R. Hawkins S. J. Neto A. I. Thompson R. C. (2008). Exploitation of rocky intertidal grazers: population status and potential impacts on community structure and functioning. Aquat. Biol. 3, 1–10. 10.3354/ab00072 Martins G. M. Thompson R. C. Neto A. I. Hawkins S. J. Jenkins S. R. (2010). Exploitation of intertidal grazers as a driver of community divergence. J. Appl. Ecol. 47:1282–1289. 10.1111/j.1365-2664.2010.01876.x Merritt T. J. Shi L. Chase M. C. Rex M. A. Etter R. J. Quattro J. M. (1998). Universal cytochrome b primers facilitate intraspecific studies in molluscan taxa. Mol. Mar. Biol. Biotechnol. 7, 7–11. 9597773 Nei M. (1987). Molecular Evolutionary Genetics. New York, NY: Columbia University Press. Nei M. Tajima F. (1981). DNA polymorphism detectable by restriction endonucleases. Genetics 97, 145–163. 6266912 Nice C. C. Shapiro A. M. (1999). Molecular and morphological divergence in the butterfly genus Lycaeides (Lepidoptera: Lycaenidae) in North America: evidence of recent speciation. J. Evol. Biol. 12, 936–950. 10.1046/j.1420-9101.1999.00111.x Nylander J. A. A. (2004). MrModeltest v2. Program Distributed by the Author. Evolutionary Biology Centre; Uppsala University. Ridgway S. A. Reid D. G. Taylor J. D. Branch G. M. Hodgson A. N. (1998). A cladistic phylogeny of the family Patellidae (Mollusca: Gastropoda). Philos. Trans. R. Soc. Lond. Ser. B. Biol. Sci. 353, 1645–1671. 10.1098/rstb.1998.0316 Rohlf F. J. (2006a). tpsUtils, Version 1.37. Geometric Morphometric Software for the PC. Copyright © 2006, F. J. Rohlf. Ecology and Evolution, SUNY at Stony Brook. Available online at: http://life.bio.sunysb.edu/ee/rohlf/software.html Rohlf F. J. (2006b). tpsDig, Version 2.10. Geometric Morphometric Software for the PC. COPYRIGHT © 2006, F. J. Rohlf. Ecology and Evolution, SUNY at Stony Brook. Available online at: http://life.bio.sunysb.edu/ee/rohlf/software.html Rohlf F. J. (2006c). tpsRelw, Version 1.44. Geometric Morphometric Software for the PC. Copyright © 2006, F. J. Rohlf. Ecology and Evolution, SUNY at Stony Brook. Available online at: http://life.bio.sunysb.edu/ee/rohlf/software.html Rohlf F. J. (2006d). tpsRegr, version 1.31. Geometric morphometric software for the PC. Copyright © 2006, F. J. Rohlf. Ecology and Evolution, SUNY at Stony Brook. Available online at: http://life.bio.sunysb.edu/ee/rohlf/software.html Rozas J. Sánchez-DelBarrio J. C. Messeguer X. Rozas R. (2003). DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19, 2496–2497. 10.1093/bioinformatics/btg35914668244 Sá-Pinto A. Branco M. Sayanda D. Alexandrino P. (2008). Patterns of colonization, evolution and gene flow in species of the genus Patella in the Macaronesian Islands. Mol. Ecol. 17, 519–532. 10.1111/j.1365-294X.2007.03563.x18179442 Slatkin M. (1987). Gene flow and the geographic structure of natural-populations. Science 236, 787–792. 10.1126/science.35761983576198 Strange J. P. Garnery L. Sheppard W. S. (2008). Morphological and molecular characterization of the Landes honey bee (Apis mellifera L.) ecotype for genetic conservation. J. Insect Conserv. 12, 527–537. 10.1007/s10841-007-9093-6 Swofford D. L. (2003). PAUP*: Phylogenetic Analysis Using Parsimony, * and other methods), Version 4.10B10. Sunderland, MA: Sinauer Associates. Tamura K. Dudley J. Nei M. Kumar S.. (2007). MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0. Mol. Biol. Evol. 24, 1596–1599. 10.1093/molbev/msm09217488738 Trussell G. C. (2000). Phenotypic clines, plasticity, and morphological trade-offs in an intertidal snail. Evolution 54, 151–166 10.1111/j.0014-3820.2000.tb00016.x10937192 Weber L. I. Hawkins S. J. (2002). Evolution of the limpet Patella candei d'Orbigny (Mollusca, Patellidae) in Atlantic archipelagos: human intervention and natural processes. Biol. J. Linn. Soc. 77, 341–353. 10.1046/j.1095-8312.2002.00102.x Weber L. I. Hawkins S. J. (2005). Patella aspera and P-ulyssiponensis: genetic evidence of speciation in the North-east Atlantic. Mar. Biol. 147, 153–162. 10.1007/s00227-004-1540-2 Whitlock M. C. (2008). Evolutionary inference from Q(ST). Mol. Ecol. 17, 1885–1896. 10.1111/j.1365-294X.2008.03712.x18363667