Congruent Molecular and Morphological Diversity of Macaronesian Limpets: Insights into eco-evolutionary Forces and Tools for Conservation

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


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 Frontiers in Marine Science | www.frontiersin.org *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.
Frontiers in Marine Science | www.frontiersin.org further morphometric analysis. In the case of shells, only nonbroken and not severely eroded specimens were considered for morphometric study.

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). 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 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).
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-platespline 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 R -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.
Frontiers in Marine Science | www.frontiersin.org FIGURE 2 | 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). 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-platespline model, based on the bending energy matrix. The congruence between molecular lineage/morphology was assessed FIGURE 3 | 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. by evaluating the degree of molecular disjunction between clades and the statistical grouping both for RW and for discriminant analysis (cross-validation).

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

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.
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. 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.
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). This morphometric analysis revealed that the lineages inhabiting the Madeira and Canary archipelagos are typified Frontiers in Marine Science | www.frontiersin.org FIGURE 7 | 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. 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 Q ST and F ST 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 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). 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 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.
(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 subspecies 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).
Frontiers in Marine Science | www.frontiersin.org