Morphometric variation in two intertidal littorinid gastropods

Here we investigate spatial variation in shell shape of Littorina saxatilis , an ovoviviparous species, and Melarhaphe neritoides , a species with planktonic eggs and larvae. Populations of both species were sampled in 6 sites located along the west coast of the Iberian Peninsula. Shell shape was studied using landmark- based morphometric methods. Landmark data was used to estimate individual size and to describe shell shape. Prior to statisti- cal analysis specimens were aligned using Generalised Procrustes Analysis (GPA). Ordinations based on GPA of landmark data and using Principal Components Analysis (PCA), revealed a clear spatial segregation of sites for both species, although this was more evident for L . saxatilis than M . neritoides . Statistical analysis revealed significant multivariate variation in shape among sites and a significant effect of allometry on shape for both species. In contrast to this last result, there was no signifi- cant linear relationship between any of the first three PC axes and size in L . saxatilis but highly significant associations between the first and third PC axes and size for M . neritoides . Spatial variation in the shape of both species was primarily related to variation in the height of the apical whorls and in the width of the aperture for L . saxatilis . Further variation in shape was related to the shape of the last whorl in L . saxatilis and shell elongation and a change in aperture shape in M . neritoides .


Introduction
Quantifying phenotypic change is a crucial first step towards understanding the evolution and ecology of species (Jensen, 1988;Rice and Mack, 1991;Cleary et al., 1995Cleary et al., , 2002Hageman and Sawyer, 2006;Adams and Collyer, 2009). Benthic invertebrates make good models to investigate spatial phenotypical variation because adults are generally sessile or sedentary and subjected to habitat-specific pressures that are relatively easy to assess. At the same time benthic invertebrates display a very large range of developmental patterns, from ovoviviparity to planktonic development lasting several months for species with teleplanic larvae (Scheltema, 1971(Scheltema, , 1986. Most often morphological traits of the shell are used to investigate phenotypic variability and plasticity given that they predictably respond to ecological factors such as wave action, desiccation risk or predator attack (e.g. Vermeij, 1973;Trussel, 2000). A common garden experiment involving two Australian species (Parsons, 1997) showed that the shell of a species with planktonic larvae, Austrocochlea constricta (Lamarck, 1822), was more plastic than the shell of a direct developer, Bembicium vittatum Philippi, 1846, indicating that wide-ranging dispersers are more physiologically flexible than limited dispersers. Moreover, B. vitatum from most source populations retained their native phenotypes, suggesting genetic control of shell shape. The different basis of shell shape variability (plastic vs genetically-determined), however, resulted in similar levels of phenotypic differentiation among native populations. These findings suggest that for wide dispersers the evolution of physiological flexibility to cope with a potential distribution across multiple adult habitats may be more advantageous than one single generalist phenotype. A different study (Hollander et al., 2006), also contrasting direct and planktonic developments involving the European species Littorina littorea (Linnaeus, 1758), with planktonic larvae, and Littorina saxatilis (Olivi, 1792), an ovoviviparous species, did not find significant differences in the degree of plasticity between the two species in a common garden experiment. The authors suggested that the scope for greater plasticity in L. littorea might be constrained by unreliable cues from wave action in wave exposed shores, or by mechanical constraints associated with disruption of shell accretion under unfavourable growth patterns. These two studies highlight that life history may not be tightly related to magnitude of phenotypic plasticity, and that contrasting mechanisms of local genetic adaptation and plasticity may ultimately produce similarly high levels of phenotypic differentiation among local populations.
An additional factor that may introduce complications in the interspecific and interpopulational comparisons of phenotypic traits, including shell shape, is allometry. Common patterns of allometric growth in gastropods include variations in the apical angle that result in doming of the shell (O'Loughlin and Aldrich, 1987), and increased thickness as an apparent antipredation adaptation (Cotton et al., 2004). Former views that allometric relationships in shell-secreting organisms, including gastropods, were genetically determined, and that spatial changes in shell morphology resulted from genetic differences and evolutionary change (Reid, 1996) are giving way to the notion that local environmental factors may lead to direct phenotypic changes of shell shape in marine gastropods with age. Experimental manipulations of the exposure of the dogwhelk Nucella lapillus (Linnaeus, 1758) to crab scent showed that the perceived risk of predation resulted in increased shell thickness, apertural tooth height, lip thickness and retractability, after correction for size differences, all of which are consistent with antipredatory defences (Palmer, 1990). In turn, the manipulation of food availability in the periwinkle Littorina littorea resulted in a more globose shell in fast growing, well fed individuals. Assuming an upper constraint in shell growth by a constant calcium carbonate deposition rate, a more globose shell would accommodate a larger body given its superior internal volume to shell material ratio (Kemp and Bertness, 1984).
In the present exploratory study, we assess to what extent the shell shape of two species of intertidal littorinid gastropods, Littorina saxatilis, an ovoviviparous species (Reid, 1996) and Melarhaphe neritoides (Linnaeus, 1758), a species with a planktonic developmental time of several weeks (Lysaght, 1941;Daguzan, 1976) differs among six sites located along the west Iberian coast, as a first step to understand the patterns of gastropod phenotypic variation in this region. We use geometric morphometric analysis to describe shape changes in snails collected from pre-determined levels of six exposed rocky shores and assess whether allometry influences shape variation in both species.

Sample collecting
Samples were collected from six populations along the Iberian coast ( Fig. 1) from October 2003 to October 2004. All sites are wave-exposed rocky cliffs where Carcinus maenas (Linnaeus, 1758), a predator with a reported impact on gastropod shell morphology (e.g. Palmer, 1990), was not detected during sampling, and for which there are no records of occurrence. Another crab species, Pachygrapsus marmoratus (Fabricius, 1787), is a very common inhabitant of rocky shores around the Iberian Peninsula and of the sampled sites (Zariquiey- Alvarez, 1968). Effects on gastropod shell morphology from predation by this crab have not been reported. Although this effect cannot be totally excluded, P. marmoratus is an omnivorous species with a relatively soft shell and flat chaela tips, which should have a much lower predatorial impact on gastropods. Each site was visited once by 2 to 3 collectors that explored a 10 to 50 m stretch of coast during approximately 15 min per collector, for each species. In order to minimize shape differences due to different ecotypes, which are a characteristic feature of Littorina saxatilis and other littorinid species, sampling was conducted within the barnacle belt for L. saxatilis and above this belt for Melarhaphe neritoides. Littorina saxatilis is a highly polymorphic species displaying a range of ecotypes that differ in size, shape, ornamentation and colour (Reid, 1996). In Galician shores (Carballo et al., 2005), two main ecotypes have been described, a ridged and banded morph found among barnacles in the upper intertidal and a small-sized, smooth and unbanded form found among mussels in the lower intertidal. In the transition zone where mussel density increases towards the typical dominance in the lower intertidal, a hybrid between the two ecotypes has been described. In order to minimize shape differences due to different ecotypes of L. saxatilis, sampling for this species was conducted exclusively within the barnacle belt, above the level of the first mussels. Melarhaphe neritoides is a typical high shore species with abundance maxima above the barnacle zone and sampling was conducted above the barnacle zone. Sampled gastropods were frozen at -20° C upon arrival at the laboratory and later transferred to 95% ethanol for 24h to facilitate the removal of soft parts.

Morphometric data acquisition
A subsample of 16 individuals from each population and species was randomly selected. A photograph of each specimen and of a scale was taken using a digital camera equipped with a 105 mm macro lens. Gastropod shells grow by spiral accretion of calcium carbonate starting at the larval protoconch and ending at the labrum of the adult teleoconch (Fretter and Graham, 1994). Because the protoconch at the apex and the apex itself are the most delicate parts of the shell and are usually eroded in larger juveniles and adults, implicating that the starting point of ontogeny is normally unknown, the determination of homologous points of the gastropod shell anatomy is often impossible (Guralnick and Kurpius, 2001). Moreover, since new morphology is continuously being added as the shell grows, older individuals do not share homologous landmarks in the new part of the shell with younger individuals. To minimize this problem, shells were mounted with the columella axis aligned with the horizontal and the shell rotated around this axis so that the aperture presented a frontal view to the camera (ventral view). Landmarks corresponding to 8 anatomical features (Fig. 2) were digitised using tpsDig software (Rohlf, 2006). These landmarks were selected because they: i) can be recognised in both species; ii) retrieve significant information related to the shape of the last whorl; iii) retrieve significant information relative to the elongation of the shell, given that erosion of the apex is not very pronounced in the two species; and iv) convey information related to the size and shape of the aperture. Landmarks 1 and 7 are Type I, 2 and 4 Type II and 3, 5, 6 and 8 Type III landmarks (Bookstein, 1991). Landmark data were used to estimate individual size, which was calculated as centroid size, i.e. the square root of the sum of the squared distances between the centroid and each landmark. In the absence of allometry, centroid size is uncorrelated with measures of shape (Tang and Pantel, 2005).

Statistical analysis of shape variation
We used Procrustes superimposition to compare shape, which minimises the sum of squared distances between pairs of landmarks on two different samples by adjusting size, rotation and translation. This squared distance is known as the Procrustes distance. Claude (2008) defines a number of functions to remove the effects of size, location (e.g. the transl() function on p 149), and orientation. The comparison of several samples (or configurations) is non-trivial because it entails a generalised procedure for the superimposition of all these samples and requires the definition of a single objective reference known as the mean shape. Shape variation is then compared in reference to this mean shape. In the present study, we use Generalised Procrustes Analysis (GPA), which determines the mean shape as the shape where the sum of pairwise squared coordinates with other rotated samples is minimised (Claude, 2008). A number of iterative algorithms have been described to obtain the best fit (Gower, 1975;Rohlf and Slice, 1990). In addition to this, full GPA entails posterior scaling between samples and mean shape to improve fit, while partial GPA maintains all samples at unit size. In the present study, we used the procGPA() function in the shapes library, which performs a full GPA. The default arguments were used, which included leaving the scale and tangentresiduals arguments as TRUE. In the results, we present Principal Components Analysis (PCA) ordinations using raw scores obtained with the procGPA() function. In addition to this, we use the shapepca() function in the shapes package to provide graphical summaries for PC's of shape. Summaries are provided whereby vectors are drawn from the mean to +1 sd along the first 3 PC axes and a Thin-plate spline (TPS) grid from the mean to +1 sd along the first 3 PC axes. TPS's are a means of interpolating how shape change affects the whole shape of an organism and mathematically express D'Arcy Thompson's deformation grids (cited in Claude, 2008). The mathematics involved in TPS analysis were imported from continuum mechanics, where they were used to study the bending of thin metal plates under physical strains (Bookstein, 1989).
We tested for significant multivariate variation in shape among different populations of both snail species using a modified Hotelling-Lawley trace statistic and F-approximation in a function taken from Claude (2008). Importantly, when working with superimposed data, the shape space dimensions are not equal to the number of variables analysed. The transformations, namely, reduce the rank number of the original data leading to four lost dimensions in 2D and seven in 3D data respectively. Thus when working with Procrustes data, one needs to adapt the normal tests. To begin with, the amount of variables in the transformation of the multivariate statistic to its F approximation needs to be replaced by the number of space dimensions. A Moore-Penrose generalised inversion also needs to be applied instead of classic algorithms due to problems in inverting nonsingular matrices as opposed to normal covariance matrices.
We analysed the relationship between shape and allometry following the procedure outlined by Claude (2008). In summary, this entails first generating a linear model between shape and allometry using the lm() function in R and subsequently using the modified Hotelling-Lawley trace statistic to test for significance. In addition to the above multivariate test, we also tested for significant variation in centroid size among populations with the lm() function in R after controlling for deviations from normality with the shapiro. test() function in R (there were no significant deviations, i.e. P > 0.05, for both species) and an effect of allometry on shape using linear regression of PC scores versus centroid size. Examination of the residuals of these regressions, however, revealed both long and short tailed distributions in addition to the presence of outliers. Although square root transformation of the response variable improved the behaviour of the residuals, the fit was still less than perfect. We therefore decided to analyse the data using robust regression (Venables and Ripley, 2002). Robust or resistant regression methods dampen the effect of outliers on model fit. There are a large number of different methods with advantages and disadvantages. In the present study we used MM-estimation as proposed by Yohai et al. (1991) and implemented in the MASS library in R. MM-estimation combines the resistance of resistant regression methods with the efficiency of M-estimation (Venables and Ripley, 2002).

Results
The first three PC's accounted for 69% (33, 20 and 16% respectively) and 62% (32, 16 and 14 respectively) of the total variation in shell shape of L. saxatilis and M. neritoides respectively. Ordinations of these axes reveal a clear spatial segregation of sites for both species, although this is more evident for L. saxatilis than M. neritoides with the degree of separation more pronounced along the first and third axes than the second axis (Fig. 3). There is, however, substantial overlap, particularly among locations in the central part of the Iberian coast.
In the case of L. saxatilis (Fig. 4), PC1 summarised variation in the height of the apical whorls (relative positions of landmarks 1, 7 and 8 along the shell axis) and in the width of the aperture (separation between landmarks 2 and 4), with a narrowing of the apical whorls associated with a widening of the aperture. PC2 was associated with the shape of the right side of the last whorl (a displacement of landmark 4 towards the bottom of the shell was accompanied by an approach of landmark 2 along a perpendicular direction) and with the curvature of the labrum (displacement of landmark 8 perpendicularly to the shell axis). PC3 was also mostly related to the shape of the last whorl, but this time on the left side (migration of landmark 3 along the shell axis). Regarding M. neritoides (Fig. 5), PC1 was also related to the height of the apical whorls (relative positions of landmarks 1, 7 and 8 along the shell axis), as was the case with L. saxatilis. PC2 was associated with overall elongation of the shell (simultaneous displacement of landmarks 3 on the left and 4 and 8 on the right in a direction normal to the shell axis). Finally, PC3 was associated with a change in aperture shape (reflected by displacements of landmark 4 along the shell axis) and with a reduction in height of the apical whorls (displacement of landmarks 5 and 6 along the profile of the shell relative to the apex).

Discussion
The present study found significant differences in shell shape among populations of Littorina saxatilis and Melarhaphe neritoides along the west coast of the Iberian Peninsula. The most conspicuous differences in shape, which were common to both species and described by PC1, were related to a variation in the height of the apical whorls relative to the rest of the shell, in such a way that a decrease in height of the apical whorls resulted in a relative widening of the aperture. This configuration of the shell is able to accommodate a larger body but with a relatively smaller foot in the snails with higher spire.
The causes for this common pattern of variation are unclear. Melarhaphe neritoides lives in the higher intertidal and supratidal zones, while Littorina saxatilis has a more extended distribution occurring from the lower supratidal zone to the lower intertidal. Ecotypic variation of the shell of littorinids in rocky shores is related to wave action, desiccation risk and predation by crabs (Johannesson et al., 1993;Newkirk and Doyle, 1975;Reid, 1993). Often, more than one ecotype can coexist in the same shore at different levels, as a consequence of the steep environmental gradient associated with the vertical position on the shore, the low vagility of individuals within microhabitats that cause assortative mating (Erlandsson et al., 1998), and the absence of a planktonic dispersal phase (Johannesson, 2003). Effects of wave exposure, desiccation and predation in littorinids and other marine snails are correlated to a large extent (Elner and Raffaelli, 1980;Palmer, 1990). Wave-exposed shores have less predatory crabs and the risk of dislodgement is the main selective pressure; this tends to produce globular shells, which reduces friction, and with a large aperture to accommodate a large foot. In protected shores there is greater predatory pressure and the snails may suffer long periods of desiccation because of reduced wave splash; this environment favours thicker shells and small apertures. Because sampling of M. neritoides and L. saxatilis was restricted, respectively, to the lower supratidal and to the barnacle belt, it is not likely that ecotypic variation associated with shore level is the cause of the differences in shell shape within species. Moreover, living in different levels, the two species are exposed to different intensities of predation, wave impact and desiccation, and should show independent responses. Although the type of shape variation accounted for by PC1 was similar, average PC1 scores of the two species were uncorrelated across sites (not shown), indicating that this component of shape was different between the two species in the same shore. This finding is also consistent with independent responses associated with shore level. There was pronounced allometry in M. neritoides with regard to PC1, and to a lesser extent PC3, that was not found in L. saxatilis (although the multivariate test revealed a significant effect of allometry in the shape of L. saxatilis), thereby indicating that different processes are controlling the variation in spire height in the two species. The results show that larger M. neritoides have an increased volume in the apical part of the shell, suggesting a larger number of whorls as a consequence of faster growth or older age. It is unclear what is driving shape variation in both species but the possibility exists that in L. saxatilis it is driven by ad-aptation to microhabitat characteristics. This species is ovoviviparous and has dispersal distances during juvenile and adult phases usually < 2 m (Erlandsson et al., 1998). Taken together, these traits underlie the extended phenotypic variability of the species and the evolution of distinct morphs at the microscale that have a genetic basis (Johannesson et al., 1995;Makinen et al., 2008, Peterson andFry, 1987). All shores sampled were exposed to the open ocean, but presented topographic variations at the microscale related to the existence of crevices, folds and general angle of the rock surface facing the waves, which were not controlled for in the present study. It is possible that changes of wave impact and water flow at scales of a few meters within the same level of the shore might lead to phenotypic variation related to the proportion of body and foot masses and, consequently, of shell proportions.
Because of the north-south orientation of the western Iberian coast some clinal variation of shell shape could be expected, related to the warmer temperatures in southwest Portugal. For instance, an increase in spire height towards the equator in high intertidal species, with corresponding reduction in aperture size and contact area of the foot with the substrate, has been observed and interpreted as favouring heat loss (Clarke, 1983;Vermeij, 1973). Also, it has been suggested that the high cost of calcification in cold waters, because of the higher solubility of calcium carbonate (Clarke, 1983), would result in thin shells with a globular shape, which are more efficient in providing a large internal volume with minimum shell material (Reid, 1996). We did not observe a significant monotonic relationship between shape and latitude in the present study but there was some suggestion of a unimodal relationship in Littorina saxatilis, particularly along PC1 with higher values in the central part of the range. Along the west coast of Iberia, yearly-averaged sea surface and near surface air temperatures range from 18.5 and 19.0 °C, in the south, to 13.0 and 15.5 °C in the north. However, because of an intensification of upwelling in central and northern Portugal relative to the rest of the coast (Peliz et al., 2005), there is a zone of cooler sea surface temperatures during spring and summer extending from the border between Portugal and Spain, in the north, to S. Martinho do Porto (Lima, 2007). The balance between these two patterns results is a zone with smaller temperature variations, lacking the cooler temperatures in winter that are found to the north and where the temperature increase during summer is less marked. Additionally, this is an area of persistent reduced salinity because of river runoff. These features are associated with a distribution gap of at least one rocky shore species located in the area (Lima, 2007). It is possible that these changes in oceanographic conditions also influence the unimodal relationship of shape and latitude, but this is an issue that would need further research.
We cautiously interpret the patterns of morphological variation identified in the present study as a consequence of the life history traits of the two species. In the ovoviviparous Littorina saxatilis, low gene flow may allow local populations to evolve through genetic drift or adaptation to local conditions. Conversely, high gene flow in the case of the plankton-disperser Melarhaphe neritoides may prevent drift and morphological change. Low levels of gene flow, and the scope for geographical differences in morphology that it allows, have been implied in most of the literature on L. saxatilis (e.g. Janson, 1982;Makinen et al., 2008;Morgan, 1999;Newkirk and Doyle, 1975). Concomitantly, for M. neritoides, high gene flow has been reported across thousands of kilometers based on allozyme markers (Johannesson, 1992). It is, however, possible that the observed spatial variation in the shell shape of L. saxatilis may also be at least partially due to non-allometric plasticity, which was previously shown to be relatively high in the species (Hollander et al., 2006). The morphological differences observed may thus represent a plastic response but unrelated to size. Direct experimental evidence is thus necessary before we can be identify to what extent the morphologic differences in L. saxatilis are due to local adaptation, allometry, or non-allometric plasticity.