Genetic and morpho-physiological differentiation in a limpet population across an intertidal gradient

Study sites and sampling

This study was carried out in the Waikiki rocky shore, (hereinafter WRS) (38°04′50S, 57°30′08W), Buenos Aires, Argentina. Average WRS tidal amplitude is 0.80 m (maximum 1.69 m) and during low tide, the beach is 10–20 m wide.

Field sampling for this study was performed during April 2010 (mid-fall) at three selected shore heights. The lower shore was defined as the minimum sea level at spring tides, which was dominated by mussel beds. The middle shore was defined as the maximum sea level at neap tides, which was predominantly covered by diverse species of native algae. The upper shore was defined as the maximum sea level at spring tides, dominated by bare rock with mussels and barnacles in cracks and crevices [29]. During sample collection, patchy micro-habitats such as crevices were avoided in order to avoid confounding habitat variation with shore height variation.

Population size structure

To determine field distribution patterns of limpets, we sampled all individuals present in ten (50 × 50 cm) quadrats randomly placed along the three shore heights. Size-frequency distributions were developed for each shore height after sorting measured limpets into 1 mm size classes. A permutational multivariate analysis of variance (PERMANOVA) [33, 34] was performed to determine differences in population size-structure among shore heights. Prior to PERMANOVA, a test of multivariate homogeneity of variances (PERMDISP) [35] was run to find differences in within-group dispersion based on the Bray–Curtis index applied to fourth-root transformed data. The significance of the factor shore heights in both analyses was tested with 999 permutations. An analysis of similarity percentages (SIMPER) [36] was performed for data recorded for each size class present at each shore height. The SIMPER analysis identified the size ranges that contribute most to the observed differences for shore height by the Bray–Curtis similarities between samples.

Morphometric analyses

To analyze shell morphology among different shore heights, a set of 20 S. lessonii adult individuals were randomly collected at each shore height. We used adult specimens (shell length > 6 mm) because smaller limpet shells are easily broken by manipulation.

The shells were photographed with a digital camera (Canon Power Shot A580) from the right side. Digital images were taken against a white illuminated background in order to maximize the contrast of shell outlines. All the images were binarized (i.e. transformed into white for the shell outline and black for the background, in pixels) so the outlines of each continuous contour (interface between the black and the white pixels) were automatically obtained and digitalized using the SHAPE software [37]. The shell shape of these limpets is rather simple with very few homologous points to be used as landmarks. Moreover, the landmarks are difficult to locate, being classified as type 2 landmarks (maximum curvature along the boundary or outline of the specimen) [38]. Shell shape variation among S. lessonii individuals from the three shore heights was therefore measured using outline analyses based on the Elliptic Fourier analysis on the outline coordinates [39]. Elliptic Fourier analysis are preferred over classical morphometric analyses in cases landmarks are difficult to determine, and this analysis have been used before in S. lessonii [25]. Elliptic Fourier coefficients were mathematically normalized in order to avoid biased outcomes resulting from different sizes, locations, rotations and starting position of shells [39]. The closed curve of each shell was broken down into 15 harmonically related ellipses. These 15 harmonics represent 99.99% of the total Fourier power spectrum [40].

A PERMANOVA test with shore height as fixed factor was used to analyze the main morphometric differences. Prior to PERMANOVA, a PERMDISP test was run to find differences in within-group dispersion based on Bray–Curtis distances applied to fourth-root transformed data. For both cases, the significance of the factor shore height was tested with 999 permutations. A posteriori Tukey pairwise comparisons were subsequently conducted. Non-metric multidimensional scaling (NMDS) was used to show graphically the morphometric dissimilarities among individuals across shore height based on Bray–Curtis dissimilarity.

Water loss regulation capacity and mortality

In order to take a representative sample of all sizes and to have replacements in case they died during transport and/or acclimatization, a total 424 individuals of all the sizes found in each shore height (numbers of individuals: upper = 140, middle = 145, lower = 139) were randomly collected. A total of 405 out of 424 individuals were used in laboratory experimentation (the 19 rest were released). To minimize stress, the collected individuals were transported to the laboratory in an icebox with rocks collected in the natural environment.

Prior to the measurements, the limpets were acclimated in aquaria with a continuous seawater flux for at least 7 days. The water flow simulated a waterfall with a seawater spray similar to what limpets are exposed to in nature. Limpets were periodically feed by daily addition of rocks with biofilm in the systems. Water temperature was the same as the ambient seawater temperature (typically between 10 and 18 °C).

Evaporative water loss (W_L) was measured according to McMahon and Britton [41] and Sokolova and Pörtner [42]. Prior to starts the experiment, limpets were removed from the aquaria and blotted with tissue paper to remove excess water from the shell surface. Shell length was measured to the nearest 1 mm to separate individuals into size groups. Individuals were weighed to the nearest 0.025 mg and placed in groups of 9 individuals (3 small: 4–7 mm, 3 medium: 8–11 mm, and 3 large: 12–15 mm) in 250 cm³ plastic bowls. These plastic bowls were incubated in a thermostatic chamber (INGELAB, model I-209D) with controlled photoperiod cycle of 8 L: 16 D. The thermostatic chamber was set with the maximum range of seawater temperature registered in the field (around 18 ± 3 °C) for the sample period [43]. The thermostatic chamber generates extra humidity that may result in higher water condensation in the plastic bowls, so silica gel was added to the experimental chambers to prevent this water condensation. After exposure periods of 12, 24 or 48 h, five plastic bowls from each shore height were removed from the chambers, and the limpets were weighed (a total 45 plastic bowls). Then they were placed in seawater, left to recover for 6 h and scored for mortality. The number of individuals used in the final analysis of evaporative water loss varied depending on the number of dead limpets. In order to estimate survival during the experiment, we identified and recorded the number of dead limpets. To estimate all the parameters of water loss in the equation, live limpets used in laboratory experiment were sacrificed (by freezing) and placed in pre-weighed aluminum pans, dried for 12 h at 75 °C, and then weighed.

Data normality (Shapiro–Wilk test) and homogeneity of variances (Cochran’s test) were tested and when necessary, data were transformed to meet statistical assumptions [44]. The slopes and elevations of the regressions of water loss W_L (dependent variable) were tested with analysis of covariance (ANCOVA) to assess the effect of shore height and time of exposure on the W_L using limpet shell length as the covariate. Two-way ANCOVA analyses can be used to compare elevations of regression lines if their slopes are not statistically different [43]. When slopes were different, we used Tukey multiple comparison tests [44] to determine which combinations of slopes differ. In these cases, we applied the Johnson–Neyman test [45] to identify the range of the covariate (i.e. shell length) where the elevations are not significantly different.

A generalized linear model (GLM) was used to evaluate whether the proportion of dead limpets (dependent variable) could be explained by shore height, time of exposure and their interaction using limpet shell length as the covariate (explanatory variables). The model was fitted using binomial distribution with logit link function [46]. When slopes were heterogeneous, interaction means comparison tests were used to determine which combinations of slopes differ for GLM [47].

Genetic ISSR analyses

To analyze the genetic diversity among shore heights, an extra set of 20 individuals of S. lessonii adult specimens (shell length > 6 mm) were randomly collected at each shore height. We used adults because small limpet shells are easily broken by manipulation and could be a contamination factor for genetic analysis. Genetic diversity was estimated using the ISSR-PCR technique. Inter simple sequence repeats (ISSR) provided a new dominant genetic marker that amplifies nuclear non-coding DNA using arbitrary primers. Primers amplify DNA fragments between inverse-oriented microsatellite loci, with oligonucleotides anchored in the microsatellites themselves.

DNA was extracted from small pieces of mantle tissue using the Chelex 100 (Biorad) method sensu Walsh et al. [48]. The primers (AG)₈Y and (CT)₈GT were used due to the high polymorphic results yielded in previous tests. The amplification reaction was performed with 20 μl final volume including 2 μl of 10X buffer with MgCl₂ (1.5 mM), 1 μl of dNTPs (2.5 mM), 4 μl of each primer (2 mM), 1 μl of template and 0.08 units/ml AmplitaqTM (Sigma), completing the remaining volume with water. The PCR reaction included an initial denaturation cycle at 94 °C (2 min); followed by 5 cycles at 94 °C (30 s), 50 °C (45 s) and 72 °C (1 min) and then another 35 cycles at 94 °C (30 s), 40 °C (45 s) and 72 °C (1 min) and a final extension at 72 °C (2 min). PCR products were run in 1.5% agarose gels, stained by Ethidium Bromide, with a molecular weight marker (1 kb). ISSRs were visualized using a UV transilluminator and analyzed by digital photography. ISSR bands with high intensity were recorded as 1, while the absence of the band was recorded as 0.

To evaluate the importance of “between-groups” (among shore heights) differentiation relative to “within-group” (for each shore height) variation of the morphometric data, a multivariate analysis of variance (MANOVA) together with Principal Coordinates Analysis of distance matrix were performed with PAST [49]. A Wilk’s Lambda test was carried out to detect if there were significant morphological differences between groups; and post hoc Hotelling pairwise comparisons (Bonferroni corrected and uncorrected) were conducted to detect significant differences using PAST.

Size-structure

Morphometric analyses

Analysis of the multivariate homogeneity of variances among S. lessonii shell shapes at different shore heights was significant (F_2,52 = 8.61, P-value < 0.01) (Fig. 1a). Subsequent pairwise test showed that the upper shore was different from the middle and lower shores (P-value < 0.01). In general terms, the variation was explained mainly by shell height (Fig. 1b), with the upper shore showing less variable shell shape than the middle and lower shores.

Water loss regulation capacity and mortality

The overall results show that the individual shell size has a positive association with the loss of water (see Fig. 2). The two-way ANCOVA showed that size had a significant effect on W_L. Larger individuals from all three shore heights lost water faster under desiccating conditions in air (Table 2). We found that the changes in water loss related to shore height was dependent on the time of exposure (Table 2). Multiple comparison of 12 h showed that the slopes and elevations among shore heights did not differ. At 24 h, comparisons showed that slopes of all combinations did not differ while the elevation for lower shore was higher than for middle and upper shores. At 48 h, comparisons showed that the slope of middle shore was lower than lower and upper shores. The Johnson–Neyman test detected no difference in the elevation of the lines depicting the relationship between W_L and shell length only in limpets of 1.5–2.1 (corresponding to 4.5–8.2 mm for non-transformed data) between lower and middle shores, and only in limpets 2.3–2.75 (corresponding to 10–16 mm for non-transformed data) between lower and upper shores (Fig. 2).

Two way Ancova	df	F
Slope (shore:size)	2	3.029*
Slope (time:size)	2	11.260***
Slope (time:shore)	4	4.034*
Slope (time:size:shore)	4	0.917 n.s.
Elevation (shore)	2	27.624***
Elevation (time)	2	80.091***

The dependent variable and covariate are expressed as natural logarithms. *p < 0.05; ***p < 0.001; n.s., not significant

The overall results show that individual shell size has a logistic negative association with the rate of mortality (see Fig. 3). The GLM showed that smaller individuals were significantly more likely to die than larger individuals (Table 3, Fig. 3). We found that the changes in the proportion of dead of S. lessonii individuals at different shore heights were dependent on the time of exposure. Multiple comparison of 12 h showed that the slopes and elevations among shore heights did not differ (Table 3). At 24 h, comparisons showed that the slope of lower shore differed from the slopes of middle and upper shores, the main difference was that lower shore decreased with size more slowly than middle and upper shores (Fig. 3). At 48 h, comparisons showed that the slopes of all combinations differed significantly, the main difference was that the lower shore decreased with size more slowly than middle shore, and middle shore more slowly than upper shore (Fig. 3).

Analysis of deviance	df	Deviance
Slope (shore:size)	2	5.983*
Slope (time:size)	2	10.435***
Slope (shore:time)	4	16.978***
Slope (shore:time:size)	4	0.514 n.s.
Elevation (shore)	2	45.001***
Elevation (time)	2	86.638***

GLM was fit under binomial distribution and logit link function. *p < 0.05; ***p < 0.001; n.s., not significant

Genetic ISSR analyses

ISSRs (AG)₈Y and (CT)₈GT showed 9 and 10 loci (different bands in the agarose gel), respectively. Total genetic diversity (H) was 0.4082 (sd = 0.027) while the diversity within groups (I) was 0.3758 (sd = 0.0072). The comparison among shore heights showed that the effective allele number, Shannon index and genetic diversity were higher for lower shore than for upper and middle shores (Table 4). Global genetic difference (G_ST) was 0.0795. The AMOVA indicated that 7% (P-value < 0.05) of the genetic variability among shore heights was explained by the population structure. The genetic differentiation PhiPT between shore heights sampled pairs was mainly ascribed to the difference between lower–middle (0.099, P-value < 0.05) and lower–upper shore (0.059, P-value < 0.05) and, to a lesser extent, between middle-upper shore (0.036, P-value < 0.05). Figure 4 shows the Principal Coordinates Analysis (PCoA) of the genetic distances among shore heights. This figure illustrates that individuals from the lower shore form a group which is slightly separated from the middle and upper shores.

Shore heights	Ne	I	H
Lower	1.747 (0.065)	0.592 (0.032)	0.409 (0.027)
Middle	1.555 (0.073)	0.501 (0.038)	0.330 (0.032)
Upper	1.675 (0.058)	0.570 (0.029)	0.388 (0.025)
Total	1.659 (0.039)	0.554 (0.02)	0.376 (0.017)

Standard deviations are shown between parentheses