Disentangling the genetic and morphological structure of Patella candei complex in Macaronesia (NE Atlantic)

Abstract The uptake of natural living resources for human consumption has triggered serious changes in the balance of ecosystems. In the archipelagos of Macaronesia (NE Atlantic), limpets have been extensively exploited probably since islands were first colonized. This has led to profound consequences in the dynamics of rocky shore communities. The Patella candei complex includes various subspecies of limpets that are ascribed to a particular archipelago and has been the focus of several taxonomic surveys without much agreement. Under a conservational perspective, we apply morphometric and genetic analyses to test subspecies boundaries in P. candei and to evaluate its current population connectivity throughout Macaronesia (Azores, Madeira, and Canaries). A highly significant genetic break between archipelagos following isolation by distance was detected (F ST = 0.369, p < .001). Contrastingly, significant genetic differentiation among islands (i.e., Azores) was absent possibly indicating ongoing gene flow via larval exchange between populations. Significant shell‐shape differences among archipelagos were also detected using both distance‐based and geometric morphometric analyses. Adaptive processes associated with niche differentiation and strong barriers to gene flow among archipelagos may be the mechanisms underlying P. candei diversification in Macaronesia. Under the very probable assumption that populations of P. candei from each archipelago are geographically and/or ecologically isolated populations, the various subspecies within the P. candei complex may be best thought of as true species using the denomination: P. candei in Selvagens, Patella gomesii in Azores, Patella ordinaria in Madeira, and Patella crenata for Canaries. This would be in agreement with stock delimitation and units of conservation of P. candei sensu latu along Macaronesia.


| INTRODUCTION
Conservation efforts applied to human-exploited and threatened species require a comprehensive knowledge about population structure and factors that shape differentiation within a species (Lande, 1988).
For instance, natural and/or human-induced barriers across a species distributional range can hasten isolation between populations by restricting the amount of individuals that can freely migrate (Barber, Palumbi, Erdmann, & Moosa, 2000;Kelly & Palumbi, 2010). Such constraints to connectivity can lead to genetically deprived populations that, in the face of intense human exploitation, are likely to succumb and disappear. To overcome isolation, many species have evolved lifehistory traits that are able to maximize dispersion and connectivity across large geographical areas (e.g., White, Fotherby, Stephens, & Hoelzel, 2011). For instance, many widely distributed marine organisms with limited adult movement avoid population differentiation by exhibiting large population sizes with high levels of fecundity and by releasing larvae that can potentially disperse in the water column for a considerable amount of time until they reach their final destination (e.g., Faria, Froufe, Tuya, Alexandrino, & Pérez-Losada, 2013). However, the extent of successful dispersion, which is a major determinant of population dynamics and structure, is a function of multiple and often interacting factors. For instance, an extensively and well-mixed larval pool does not necessarily lead to widespread connectivity and lack of population genetic structure over large spatial scales (e.g., Keever et al., 2009). The complex interplay between physical processes (e.g., coastal topography, stratified water columns, tidal forces, wind, buoyancy, surface waves, and turbulence) and life-history traits (e.g., time of spawning, larval behavior, growth and survival rates, pelagic larval duration), often interacting at fine to mesoscales, can result in a broad range of dispersal and metapopulation connectivity patterns (see review in Cowen & Sponaugle, 2009). Similarly, historical events such as past glaciations and changes in sea level can determine the contemporary distribution of populations (e.g., Portnoy et al., 2014).
In the last few decades, genetic methods have become a tool of excellence in investigating population differentiation estimates of a given species across its distributional range. Particularly, highresolution nuclear markers such as microsatellites have been widely applied to conservation genetics (Guichoux et al., 2011) and to the identification of populations requiring prioritizing protective measures (e.g., Sandoval-Castillo & Beheregaray, 2015).
The Macaronesia comprises five NE Atlantic archipelagos: Azores, Madeira, Selvagens, Canaries, and Cape Verde. The region is defined as a biogeographical entity based on the existence of many shared elements in the flora and fauna among archipelagos. All are of volcanic origin and have distinct but fairly recent geological times of origin (Ávila et al., 2016). Patellid limpets inhabiting these archipelagos are considered a valuable resource and have been intensively exploited presumably since islands were first colonized (Côrte-Real, Hawkins, & Thorpe, 1996;Hawkins, Côrte-Real, Pannacciulli, Weber, & Bishop, 2000;Santos, Hawkins, Monteiro, Alves, & Isidro, 1995). In most islands, heavy exploitation has led to dramatic decreases in limpet abundances with current populations showing clear signs of over-exploitation (Martins, Jenkins, Hawkins, Neto, & Thompson, 2008;Martins et al. 2017). Because limpet's grazing activity acts as a key process in shaping the structure and functioning of rocky shore communities (Hawkins & Hartnoll, 1983), the chronic removal of limpets has led to an upward extension of turf-forming algae (see Boaventura et al., 2002;. In fact, under reduced numbers of limpets, algal spores can opportunistically grow to a size that allows them to escape grazing and thus form mature algal patches that are able to persist through time altering the community dynamics and energy flow (Coleman et al., 2006;. The intertidal limpet P. candei, which is exclusive to Macaronesia, occurs on rocky shores from the mid intertidal down to 5 m depth across all archipelagos except in Cape Verde where it is absent (Christiaens, 1973). They are broadcast spawners with fertilization occurring in the water column. According to Martins, Santos, and Hawkins (1987), P. candei is a gonochoric species that spawns throughout the year, without synchronized resting periods. While very few studies have sentenced the pelagic larvae duration (PLD) in species of the genus Patella (Dodd, 1957;Ribeiro, 2008) with gross estimates ranging from 2 to 32 days depending on the species, temperature, and settlement cues, there is still no available information about the PLD of P. candei. Individuals of P. candei have a suboval to stellate shell shape and an orange to grayish foot with a thin darker border. Morphological plasticity associated with specific micro-habitat conditions (i.e., substrate complexity) and environmental variation (i.e., wave exposure) is known to occur in this species (Hawkins, Côrte-Real, Martins, Santos, & Martins, 1990). For instance, two distinct habitat-related morphs of P. candei are referred for the Azores: the "fly limpet" and the "smooth limpet" (Hawkins et al., 1990). Moreover, the morphological variation associated to each archipelago has led Christiaens (1973) to describe four distinct subspecies within P. candei complex: P. candei gomesii in Azores, P. candei ordinaria in Madeira, P. candei candei in Selvagens, and P. candei crenata in the Canaries. This classification is not entirely supported by subsequent studies. For instance, Côrte-Real et al. (1996) found no differences in radular morphology and soft-body parts among archipelagos but showed that shell shape and allozyme characters from P. candei in Azores were clearly distinguished from P. candei in the Madeira and the Canaries. Weber and Hawkins (2002) also showed that P. candei shell shape could be distinguished among archipelagos and that allozyme retrieved two well-differentiated groups: P. candei from Azores and Selvagens, and P. candei from Madeira and Canaries. More recent research, using mtDNA, showed that samples of P. candei from the Azores, Madeira, and Desertas (located about 25 km southeast of Madeira Island) form a well-supported group, while individuals sampled in Selvagens and Canaries always grouped together in a different clade (Sá-Pinto, Branco, Harris, & Alexandrino, 2005;Sá-Pinto, Branco, Sayanda, & Alexandrino, 2008). The taxonomic status of Macaronesian limpets thus remains unclear with much controversy as to whether P. candei from the different archipelagos should be given a specific status.
In this study, we use morphometric and molecular genetic methods to assess the existence of distinct groups of P. candei across Macaronesia, testing the subspecies boundaries within the species.
At multiple spatial scales, we evaluate the degree of contemporary connectivity among populations and hypothesize that populations geographically closer to each other are likely more related and connected via larval dispersion. Besides the assessment of genetic diversity and structure of P. candei populations across archipelagos, and given the importance of defining conservation units in fisheries planning , we provide discussion and guidance about protective measures of such threatened marine resource, highlighting the importance of considering levels of genetic diversity in populations as well as their uniqueness.

| Distance-based analyses
Individuals of P. candei were collected across the Macaronesia archipelagos of Azores, Madeira, and Canarias, in 12 different islands (n = 917) during the summer of 2011 ( Figure 1). In all islands, individuals were collected on intertidal platforms. The soft tissue of all individuals was carefully removed for genetic analyses, and the shells were marked individually. Shell morphometry was examined using the procedure described in detail by Cabral (2007). In summary, five distances were measured on each shell: shell length (SL), shell maximum width (SW), shell width at apex (SWA), distance from apex to anterior tip (SAA), and shell height (SH). These distances were then used to calculate base ellipticity (BE), base eccentricity (BEC), conicity (CO), and cone eccentricity (CE) for each individual shell (see Cabral, 2007 for further details). Shells were measured using a Vernier caliper with a precision of 0.1 mm. Individuals with signs of shell damage or with clear home fitting deformation were discarded from the analyses. The four variables were inspected for collinearity and removed from the analyses when collinearity was high (r > .3).
Examination of the spatial variation in shell morphometry among archipelagos was analyzed using a two-way fully hierarchical analysis of variance (ANOVA) with the following factors: archipelago (fixed) and island (random and nested within archipelago). Analyses were performed using PERMANOVA+ (Anderson, Gorley, & Clarke, 2008) based on Euclidean distances and using 999 permutations.
PERMANOVA is a permutational-based test and produces results analogous to ANOVA when based on Euclidean distances (Anderson, 2001). This test was preferred over classical ANOVA because it allows the use on unbalanced designs (e.g., different number of islands per archipelago). When needed, SL was used as a covariate in order to adjust for differences in morphometry with animal size. Prior to analysis, PERMDISP was used to test for homogeneity of multivariate dispersions. All analyses were performed on untransformed data.
Principal component analysis (PCA) was performed to determine the linear combination of morphometric descriptors that account for the variation in the data. Shell-shape variation between P. candei morphotypes ("fly" and "smooth") from Azores was also accessed by means of ANOVA. Given the difficulties in sorting all individuals into the "fly" and "smooth" categories, only distinctive shells of each morphotype were used in the analysis (n = 346). Note, however, that this criteria may upward bias the degree of differentiation between morphotypes.

| Geometric morphometric analyses
A total of 83 P. candei collected across archipelagos were used for landmark-based geometric morphometric analysis (Bookstein, 1991; Figure 1). For imaging, shells were positioned in a light background and digital high-resolution images of the dorsal surface were captured F I G U R E 1 Map of sampling locations for Patella candei collected from the Macaronesian archipelagos of Azores, Madeira, and Canaries (NE Atlantic). N columns indicate the number of individuals used for distance-based morphometrics, geometric morphometrics, and genetic analyses, respectively using a CANON EOS 600D camera mounted on a tripod to maintain the distance for all samples and to ensure that the lens was parallel to the surface examined. The anterior-posterior axis of each specimen was identified using the scars from each individual body left on the ventral surface of the shell. A fan was used to position each specimen along such axis and ensuring that the apex of each shell coincided with the vertical line of the fan (Fig. S1). A tps file of each archipelago specimens was created using tpsUtil, and tpsDig2 (Rohlf, 2015) was used to place a total of 37 landmarks on the shell apex and on the intersection of the fan and shell of each sample specimen image. Except for the shell apex landmark, all other points (at the shell border) do not necessarily represent homologous landmarks from a development point of view, but can be used to decompose objectively the shell shape of limpets. These points are referred to as semi-landmarks and can be used to capture information about curvature (Gunz & Mitteroecker, 2013). Specimens were then aligned using a Generalized Procrustes Analysis (GPA) (Rohlf & Slice, 1990) to remove all the differences due to translation, rotation, and scale (Bookstein, 1991). In this process, semi-landmarks are allowed to slide along their tangent directions so as to minimize the bending energy between each specimen and the reference form. The resulting aligned Procrustes coordinates represent the shape of each specimen.
For more details on geometric morphometric methodologies using landmarks, see Bookstein, 1991;Zelditch, Swiderski, & Sheets, 2012;Adams, Rohlf, & Slice, 2013;. Centroid size (CS), given by the square root of the sum of squared distances of a set of landmarks from their centroid, was also calculated (Rohlf & Slice, 1990). A Procrustes permutation analysis of variance (Procrustes ANOVA) performed with a residual randomization permutation procedure (Adams, Collyer, Kaliontzopoulou, & Sherratt, 2016;Collyer, Sekora, & Adams, 2015) was used to determine patterns of shell shape variation between archipelagos. The aim is to test groups (archipelagos) considering the influence of phenotypic change associated with body size (i.e., body size was calculated from landmark configurations as centroid size). Actually, the Procrustes ANOVA between shape and size, using CS is a way to assess interspecific allometry (Villegas, Feliciangeli, & Dujardin, 2002). Because allometry was detected among archipelagos (see section 3), the full dataset was divided in two: one with small limpets (SMALL) and one with big limpets (BIG).
Principal components analysis (PCA) was used to provide a graphical depiction of patterns of shape variation across the two datasets. Thin plate splines were used to provide a visual representation of the shape changes between each group mean and overall consensus configuration. All analyses and graphical representations were performed in R (R Core Team 2014) using the packages GEOMORPH (Adams & Otarola-Castillo, 2013). Geometric morphometric analysis was also used to determine shell shape variation between the "fly" and "smooth" P. candei morphotypes from Azores (n = 23 and n = 34, respectively).

| Sampling and laboratory protocols
A total of 560 individuals of Patella candei from the three archipelagos were used for genetic analysis (Figure 1). Upon collection, limpets were preserved in 96% ethanol and frozen for later processing. At the laboratory, samples were subject to DNA extraction from the foot muscle tissue using the E.Z.N.A. Mollusc DNA extraction kit and following the manufacturer' instructions. The quality and quantity of DNA extractions were assessed using a Nanodrop spectrophotometer (Thermo Scientific). All individuals were genotyped at 12 microsatellite loci using the primer pairs and following the amplification protocol described in Faria et al. (2016). Briefly, microsatellites were amplified in three distinct multiplex PCRs (PcaMix1: loci CAN18, CAN25, CAN27, CAN53; PcaMix2: loci CAN9, CAN26, CAN32, CAN40; PcaMix3: loci CAN23, CAN33, CAN56, CAN60) on 10 μl reactions containing ~30 ng DNA template, 1 × Qiagen TM Multiplex PCR Kit, 0.5-1.2 μmol/L of each primer and ddH 2 O. Genotyping was performed on an ABI 3730 (Applied Biosystems) automated DNA sequencer using an internal size standard (GeneScan TM 500Liz ® , Applied Biosystems) for accurate sizing and GENEMAPPER TM v.4.2 (Applied Biosystems) was used for allele calling.

| Genetic differentiation and population structure
Pairwise F ST estimates among populations were calculated using FSTAT v.2.9.3 (Goudet, 1995), and departures of F ST from the null hypothesis of panmixia were evaluated via a permutation test (1,000 iterations). The effect of null alleles in F ST estimates was assessed by comparing F ST before and after correction for null alleles using the excluding null alleles (ENA) method implemented in FREENA. Genetic differentiation between populations was also determined using the D est estimator (Jost, 2008) implemented in the R package DEMETICS v.0.8.4 (Gerlach, Jueterbock, Kraemer, Deppermann, & Harmand, 2010), and p-values were estimated by bootstrap analysis (1,000 replicates). For all analyses involving multiple tests, significance levels were adjusted by the FDR method.
The model-based approach implemented in STRUCTURE v.2.3.3 (Pritchard, Stephens, & Donnelly, 2000) was used to identify the most likely number of populations (K) and assign individuals to genetic clusters.
Assignment is conducted in ways that minimize deviations from Hardy-Weinberg and linkage equilibrium within each cluster. No particular population structure was assumed a priori (LOCPRIOR = 0), and ten independent runs were carried out for each value of K (1-11). Length of the burn-in period was set to 1 × 10 5 followed by 5 × 10 5 Markov chain Monte Carlo (MCMC) iterations. Correlated allele frequencies and admixed populations were assumed. Modifications in such parameters produced consistency and did not change the final results. Selection of the most likely number of genetic clusters (K) was based on checking the posterior probability of the data for a given K (Pritchard et al., 2000) and also by looking at the second-order rate of change in probability between successive K values as described in Evanno, Regnaut, and Goudet (2005) and implemented in STRUCTURE HARVESTER (Earl & vonHoldt, 2012).
In systems with hierarchical population structure, STRUCTURE typically best resolves the highest level of population subdivision (Evanno et al., 2005). Thus, in order to resolve lower levels of subdivision, structure analyses were also conducted separately for each cluster identified.
Therefore, two additional STRUCTURE analyses using the same settings were used to identify potential within-cluster structure. The best K was determined as previously described.
A discriminant analysis of principal components (DAPC) was also performed to identify and describe clusters of genetically related individuals (Jombart, Devillard, & Balloux, 2010). DAPC has been shown to perform generally better than STRUCTURE at characterizing population subdivision (Jombart et al., 2010). DAPC is a multivariate analysis that integrates principal component analysis (PCA) together with discriminant analysis to summarize genetic differentiation between groups. DAPC is free of assumptions about Hardy-Weinberg equilibrium or linkage disequilibrium and provides graphical representation of divergence among populations. DAPC was performed with and without using prior group information using the R package ADEGENET (Jombart, 2008). All STRUCTURE and DAPC analyses were conducted upon removal of nonamplifying loci.
Tests for genetic differentiation among archipelagos were also conducted using analysis of molecular variation (AMOVA) in ARLEQUIN v.3.5.1.3 (Excoffier & Lischer, 2010). Genetic variation among archipelagos (F CT ), among populations within archipelagos (F SC ) and within populations (F ST ) was assessed, and significance of F-statistics was tested using 10,000 permutations. Estimates of genetic differentiation were also determined among the "fly" and "smooth" morphotypes of P. candei from Azores.

| Isolation by distance and gene flow
To test for isolation by distance (Wright, 1943) was regressed onto the natural log of geographic distances (GD; Rousset, 1997). Regression with GD was also performed with the differentiation estimator D est matrix. Regression analyses were performed in R and tested for significance with a Mantel permutation procedure. Moreover, given the heterogeneous nature of samples, the Monmonier's maximum difference algorithm implemented in BARRIER v.2.2 was used to highlight geographical features associated with genetic discontinuities among populations (see Manni et al. 2004 for method details). Analyses were conducted using pairwise F ST values, and statistical confidence for each identified barrier was evaluated using 100 bootstrap replicates that were simulated using the package diveRsity in R. Analyses were also conducted separately for each amplifying microsatellite locus.
Analyses were only conducted among archipelagos due to the lower accuracy of BAYESASS when migration rates are high and genetic differentiation is low (see section 3) (Faubet, Waples, & Gaggiotti, 2007;Meirmans, 2014). The program was run for 3 × 10 6 MCMC iterations with sampling at every 1,000 iterations, of which 10 6 iterations were discarded as burn-in. Delta values for allele frequency, migration rate, and inbreeding were adjusted so that the accepted numbers of changes were 40%-60% of the total number of iterations. Ten MCMC runs with different initial seeds were carried out in order to maximize convergence and mixing. The Bayesian deviance was used as an optimality criterion to find the run with the best fit (Faubet et al., 2007). Deviance was calculated from the trace file using the R-script provided by Meirmans (2014).
Contemporary gene flow was also estimated using the F-model on BIMr program (Faubet & Gaggiotti, 2008). This software can estimate migration rates and detect migrants (within the last generation) at a lower level of population differentiation compared to BAYESASS (Faubet & Gaggiotti, 2008). In addition, BIMr can identify the environmental factors that are more likely to explain the observed patterns using a generalized linear model. The method employs a Bayesian approach and Markov chain Monte Carlo (MCMC) techniques to make inferences of recent gene flows in subdivided populations (Faubet & Gaggiotti, 2008). Preliminary trials included all populations but because population-specific F ST values below 0.01 can be problematic for parameter estimation, analyses were performed on samples grouped according to the Bayesian clustering analyses results. Also, analyses were conducted with and without removing loci that failed to amplify in some populations and/or exhibited null alleles. Often considered one of the main factors in determining gene flow in many species, the geographic distance between samples was included as the environmental variable. A total of 20 independent replicate runs were performed. Each MCMC was run for a total of 3.53 × 10 6 iterations, which included 30 short pilot runs of 1,000 iterations each in an effort to obtain acceptance rates between 25% and 45%. The next 15 × 10 5 iterations were discarded as burn-in, and a total of 20,000 samples were collected from each of the 20 replicates using a thinning interval of 100 iterations, using default settings. The posterior probabilities were evaluated for the run with the lowest Bayesian deviance (given by the assignment component of the total deviance: D assign ) (Faubet & Gaggiotti, 2008;Faubet et al., 2007). The mean, mode (point estimate), and 95% highest posterior density intervals (HPDI) for migration rates were recorded.

| Distance-based morphometrics
Shell samples of P. candei ranged in size (SL) between 1.35 and 6.35 cm, with a mean size of 3.25 ± 0.03 cm (mean ± SE). Although variable across islands (ANOVA p < .001), mean SL did not differ among archipelagos (ANOVA p > .05, Fig. S2, Table S1) but was significantly correlated with the remaining distance measures (SW, SWA, SAA, and SH) and also with BE (Table S2). Hence, only SL was used for analyses and considered as a covariate in the analysis of spatial variation of BE. Among the morphometric descriptors, conicity (CO) was positively correlated with BEC and CE and was therefore selected for analysis together with BE. Significant variation in the shape, as given by analyses of the descriptors CO and BE, was found at the scale of islands and among archipelagos (ANOVA p < .05, Tables S3 and S4).
Variation in shell BE among archipelagos, although affected by size (SL), was relatively higher than variation among islands (Figure 2).
Similarly, for conicity, the largest proportion of the variability was found at the scale of archipelagos (Figure 2). The first principal component (PC1) described 62.2% of the total variation, with the remaining variation (%) being accounted by the second principal component (PC2). The most important variable integrated by the first and second components was conicity and BE, respectively. The PCA showed that shell conicity can better distinguished shells from different archipelagos, whereas shell BE was mostly associated with differences within archipelagos ( Figure 3). Significant shell shape variation was also detected between the two P. candei morphotypes in Azores (ANOVA, p < .01). Differences were only detected for shell conicity, with "fly limpets" being more conical then "smooth limpets."

| Geometric morphometrics
The Procrustes ANOVA analysis on the full shell shape dataset revealed a significant interaction between archipelago and size, indicating an allometric growth in P. candei (Table 1). Moreover, the null hypothesis for common allometries (parallel slopes) among archipelagos was rejected (F = 5.364, p < .001). The amount of shape change per unit of size change differed among archipelagos and was greater in Azores (indicated by the lengths of slope vectors) (Table 1; Fig. S3).
Yet, shape trajectories and the way shapes change were only significantly different between Canaries and Madeira (indicated by the angles between slope vectors) ( Table 1; Fig. S3). This corresponds to contrasting local deformations in particular parts of the landmark configuration associated with size change in these two archipelagos (e.g., note the changes in the shell apex landmark in Madeira and Canaries; Fig. S4). Significant differences in shell shape unrelated to size were detected among archipelagos for the two subsets (SMALL and BIG; Table S5). For both datasets, pairwise comparisons showed that Azorean and Canaries shells could not be distinguished (Table S5). The first two principal components of the Procrustes shape variables for each dataset accounted for 52 and 53% of the total sample variation, respectively (Figure 4). A generalized overlapping in the scatter of data was found, mostly between Azores and Canaries samples. Intraspecific variance was greatest in smaller shells from Madeira and Canaries, with individuals from these archipelagos occupying a much wider range of shape space than samples from Azores. Deformation grids for both SMALL and BIG datasets indicate that shell goes from a clear round shape to a more ridged and pointy look-alike shape along CV1 (left to right) (Figure 4). On the same direction, the shell apex gets closer to the anterior margin of the shell. Similarly, the anterior end of the shell gets narrower along such axis. These shape changes are mostly associated and illustrate shell shape differences between Azores/Canaries and the Madeira samples. Whereas Azorean shells are oval with a smoother margin, the Canaries samples exhibit some ridges along their shell border. The pentagon look-alike shape of P. candei in Madeira stands out from the remaining archipelagos (Fig. S5). As for shell shape variation between P. candei morphotypes ("fly" and "smooth") from Azores, the Procrustes ANOVA analysis revealed a significant interaction between morphotype and size, which is indicative of allometric F I G U R E 2 Components of variability for shell conicity (CO) and shell base ellipticity (BE) in Patella candei across archipelagos; SL stands for shell length. Significance: *p < .05, ***p < .001 (see Tables  S3 and S4 for ANOVA terms) growth. The null hypothesis for common allometries (parallel slopes) among morphotypes was rejected (F = 7.459, p < .01), and differences were detected in the amount of shape change per unit of covariate change (size) between morphotypes; shape change per size unit is higher in "fly" limpets (Table S6). Overall shape of P. candei in Azores is oval for both morphotypes but the shell apex in "fly" limpets tends to get closer to the anterior margin of the shell (Fig. S6).

| Genetic analysis
A total of 138 alleles were observed across the 12 loci examined, ranging from six in CAN26 and CAN56 to 22 alleles in CAN18 (Table S7).
Five loci failed to amplify in individuals from Madeira and Canaries: Loci CAN9, CAN26, and CAN33 did not amplify at all in these two populations and two additional loci failed to amplify in more than 30% of both samples (Table S8). Since microsatellite markers were developed using P. candei from the Azores (see Faria et al., 2016), such amplification failure suggests a high genetic differentiation between Azores and the remaining archipelagos. Multilocus mean allelic richness with rarefaction was similar across populations and ranged from 3.6 (GCA) and 4.7 (FAI). Mean number of private alleles was greater in the MAD population (Ap[30] = 1.14). Observed heterozygosity frequencies (H O ) were relatively low and ranged from 0.157 to 0.364, while expected (H E ) heterozygosity frequencies ranged from 0.304 to 0.484 (Table S7). No significant linkage disequilibrium was detected for any pairs of loci. Except for loci CAN23 and CAN53, all other loci deviated from Hardy-Weinberg equilibrium (HWE). Significant locus-specific inbreeding coefficients (F IS ) ranged from 0.073 to 0.746 denoting a heterozygous deficit in such loci. Overall, significant F IS values were often ascribed with the presence of null alleles. To check for any bias in the results, loci with a presence of null alleles >10% were removed from the analysis. Such removal did not affect genetic differentiation results (data not shown) and unless stated otherwise, all loci were included in subsequent genetic analyses. In fact, the influence of null alleles has been shown to be marginal when compared to other factors such as number of loci and strength of population differentiation (Carlsson, 2008). Pairwise comparisons of F ST and D est indicated high and significant genetic differences among archipelagos but not within archipelago (i.e., between islands in Azores) ( Table 2). Both indices were highly correlated (Pearson's correlation 0.99, p < .001) (Fig.   S7). Furthermore, F ST values before and after correction for null alleles using the ENA method did not differ considerably (Table S9)    DAPC were consistent with those obtained using AMOVA, which de- The results of the migration rates estimated in BAYESASS suggest a consistent restriction in contemporary gene flow between archipelagos (Table 3). Yet, from the individual assignments, two individuals sampled in Madeira had posterior probability estimates (0.13 and 0.44, respectively) as second-generation migrants from Canaries. One individual sampled in the Canaries had a probability of 0.24 to be a firstgeneration migrant from Madeira (Table S10). Estimates of current gene flow rates from BIMr for populations of P. candei were consistent across all preliminary trials and runs and showed mean values very close to 0 (

| Shell shape variation between archipelagos
There have been few attempts to describe and distinguish limpet species and/or morphotypes using shell morphometry (e.g., Cabral, 2007;Denny, 2000). The difficulties of such methods rely on the fact that limpet shells have a suboval shape without any clear external homologous landmarks (except for the shell apex) or readily identifiable morphological features. Furthermore, due to the strong influence of certain environmental factors (i.e., wave exposure, substrate complexity, predation) many organisms, including limpets, exhibit high morphological plasticity (Branch & Marsh, 1978;Guerra-Varela et al., 2009;Harley, Denny, Mach, & Miller, 2009;Lowell, 1986;Sokolova & Berger, 2000).  Douglas, 2016;Ruane, 2015). In such sense, whereas distance-based methods are considered rather more simplistic in detecting shape differences between samples, the geometric analysis allows a more detailed recognition and assessment of how such shape varies with shell growth. Despite intrinsic different, generally, both methods allowed to distinguish P. candei morphotypes among and within archipelagos.
The differences found between both methods are likely associated to (1) the fact that distance-based data contain relatively little information about shape because many of the measurements overlap and/or are correlated to each other, and shape can only be derived from ratios among particular measurements and (2) the inclusion of shell height to the distance-based method provides a third dimension, which is absent in geometric morphometrics that only considers a two-dimensional reduction of shape in current analyses. Depicting the results, distance-based methods show that samples from Azores exhibit a more conical and elliptical shell shape compared to samples from the Canaries and Madeira, with the later descriptor being influenced by size. In fact, shell shape allometries derived from geometric morphometrics differ among archipelagos, and globally, differences are mostly found between Madeira samples and the remaining archipelagos; both small and large limpets from Azores and Canaries are more similar between them than with samples from Madeira. Whereas genetic data suggest a closer relationship between P. candei from  Trussell, 2000), can also determine geographic variation. Within the Azorean archipelago, such phenotypic plasticity is evident in the two well-recognized habitat morphs: The "fly limpet" which is highly conical and commonly found upper on the shore, mostly on more rugose surfaces; and the "smooth limpet," which is more flattened and associated to surfaces highly exposed to hydrodynamic forces (Hawkins et al., 1990). Under particular circumstances, such phenotypic plasticity can set the baseline for sympatric speciation and evolutionary divergence of habitat morphs (see Agrawal, 2001). In fact, environmental stress gradients in coastal intertidal habitats related to heat, desiccation, salinity, and wave action can provide the adequate setting for adaptive processes in patellids (Branch, 1981). If reproductive isolation is enforced by the ecological characteristics of each habitat, then biological separated species can be revealed. A good example comes F I G U R E 6 Regression between genetic distances F ST /(1 − F ST ) (top plot) and D est (bottom plot), with natural log geographical distances. Whereas a strong signal of IBD is observed among archipelagos, basal pairwise points in both graphs indicate that genetic differentiation between geographically separated islands within archipelago is absent (i.e., Azores) from the diversification in Nacella limpets in the Magellanic Province

| Population genetic structure and contemporary connectivity
Genetic differentiation estimates among archipelagos revealed a highly structured pattern in the P. candei complex across  (95% HPDI) and Canaries, which are about 400 km apart, also show significant genetic differentiation and limited contemporary connectivity. In fact, migratory events between Madeira and Canaries are unlikely, if not entirely absent, despite the fact that these archipelagos are geographically closer to each other. Only three individuals showed a very slightly probability of being migrants between these two archipelagos. Selvagens islands, which stand at approximately two-thirds of the way between Madeira and the Canary Islands, and could act as a putative stepping stone for gene flow requires further examination. As for Azores, despite the wide geographical distribution of its islands across (~600 km), the minimum and maximum distances among any pair of adjacent islands that pelagic larvae must travel among islands are approximately 32 and 220 km, respectively. Such distances do not seem to offer an obstacle and allow populations' gene pool across all islands to be homogenized via larval transport.
Although genetic differentiation among archipelagos seems to be highly correlated with geographic distance, the unbalanced nature of sampling, the fact that connectivity within archipelago (i.e., Azores) does not follow IBD, and the results provided by BARRIER (see increase almost exponentially as they move away from the coast into offshore waters (Cowen, Lwiza, Sponaugle, Paris, & Olson, 2000), but the PLD of P. candei may also be shorter than expected, or of a narrower range of what is generally referred for other patellids (Ribeiro, 2008). Under such scenario, because larvae are less likely to travel longer distances, populations farther away from each other would be more genetically distinct.
The failure of some microsatellites, which were isolated from the genome of P. candei samples from Azores, in amplifying individuals from Madeira and Canaries may further suggest substantial genetic break among archipelagos. The reduced marker polymorphism in southern samples may reflect pronounced sequence differences among subspecies due to their genetic divergence, thus entailing an upward ascertainment bias as markers were developed from Azorean morphotypes. This may explain the contrasting results of this study and those provided by Sá-Pinto et al. (2005. In their study, samples from Azores and Madeira grouped together in a wellsupported clade, leaving Canaries more distant related. According to the same authors, the scenario of a single colonization event for each archipelago associated with the absence of historical gene flow between them is the most likely. The direction and timing of such colonization events is still unclear, and methods such as the approximate Bayesian computation (Beaumont, Zhang, & Balding, 2002) (Hoskins, Higgie, McDonald, & Moritz, 2005). In fact, allopatric speciation under restrictive gene flow is believed to be one of the most common modes of speciation in nature (Schluter, 2009). Yet, providing that reproductive isolation is complete, secondary events of colo-

| Conservation of limpet population in Macaronesia
The aim of conservation is not simply to safeguard species from going extinct, but also to guarantee that morphological and genetic variation in natural populations is preserved. Efforts toward ensuring the conservation of limpet populations in Macaronesia are highly recommended, especially considering that their low genetic diversity and lack of gene exchange between archipelagos suggest they may be highly vulnerable. Taking into account its endemic nature and the negative impact of over-exploitation in coastal communities, the risk of complete extinction of P. candei in Macaronesia is therefore conceivable.
Yet, because our study failed to detect the occurrence of population bottlenecks, which would be expected under the known demographic decline of P. candei in Macaronesia, the effective population sizes must be still large enough to prevent critical losses to genetic variability (see Pujolar et al., 2011). However, when comparing to unexploited populations of patellids elsewhere (e.g., Perez et al., 2007;Ribeiro, Branco, Hawkins, & Santos, 2010), genetic diversity in the P. candei complex is fairly low. Despite the challenges associated with such comparisons, especially because different microsatellite loci may generate different levels of variation, the reduced genetic diversity observed can be a consequence of high levels of population inbreeding caused by intensive exploitation. As harvesting is mainly aimed at larger individuals (Martins et al., 2008), the more fertile individuals with higher reproductive outputs are likely less abundant than expected. This may lead to severe evolutionary and ecological consequences for the biology, life-history traits, and survival of the species (Fenberg & Roy, 2008).
Our results support the view that populations from each archipelago should be managed for conservation as distinct units. Within the P. candei complex, P. candei candei is considered in danger of extinction under the Spanish Catalogue of Endangered Species and is now mainly restricted to one single island: Fuerteventura (Canary Islands). This is likely a consequence of overexploitation (Núñez, Brito, Riera, Docoito, & Monterroso, 2003), but selective and evolutionary-related processes may have also been involved (González-Lorenzo et al., 2016). Recently, a forced ban to its capture and a conservational plan has been put in place by regional authorities. Aside from Fuerteventura, P. candei candei occurs more abundantly in the Selvagens, uninhabitable islands that are protected under the Portuguese status of Nature Reserve. As for the remaining limpets, each Regional authority has established protective measures and rules for their harvesting. Legislation does not differ much among regions, with the establishment of limpet no-take areas, minimum legal catch sizes, and seasonal fishing closures. Unfortunately, these actions have been largely ineffective in protecting such resource (Diogo, Pereira, & Schmiing, 2016;López et al., 2012;Riera et al., 2016), mostly because of illegal harvesting and lack of or insufficient enforcement. Protective measures need to be adjusted to the particular life-history traits of P. candei in each archipelago (e.g., temporal variation in reproduction, recruitment, population dynamics). Furthermore, both environmental awareness programs to general population and coastal enforcement by local authorities should be stimulated. Above all, we suggest the establishment of fully enforced closed zones to limit the access to rocky shore poachers and allow limpet populations to grow in numbers and individual sizes. These off-limit areas should be regularly surveyed and take into consideration population connectivity patterns within archipelagos. Our data suggest that, at least for Azores, such areas should distance no more than ~200 km, to allow proper gene exchange between populations. Additional sampling throughout all Macaronesia islands can further sharpen our knowledge about connectivity patterns in P. candei and help better defining the establishment of such areas within each archipelago.
Under the very probable assumption that P. candei from each archipelago forms a geographically and/or ecologically isolated population, the various subspecies within the P. candei complex as previously proposed by Christiaens (1973) may be best thought as true species using the denomination: Patella candei in Selvagens, Patella gomesii in Azores, Patella ordinaria in Madeira, and Patella crenata for Canaries. Whether this would facilitate current taxonomic misinterpretations and conservation needs, further research is still required, especially in diagnosing intrinsic reproductive isolation (Frankham et al., 2012). Ultimately, this elevation of subspecies to species level would be in agreement with stock delimitation and units of conservation  with potential benefits for management plans aimed at the preservation of limpets stocks across the Macaronesia region.

ACKNOWLEDGMENTS
We sincerely thank two anonymous reviewers for their useful comments that improved significantly this manuscript. We thank

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
JF and GMM conceived the study and all co-authors equally contributed to its final design. JF performed all the analyses. AP and PP assisted JF with the validation of genotypes and genetic analyses. GMM assisted JF with multivariate analysis. AIN coordinated sampling. JF and PR were involved in sample collection. JF interpreted the results with input from GMM, SJH, PP, and AIN. JF wrote the manuscript with input from all co-authors.