Marine Biology (2002) 140: 1215–1225 DOI 10.1007/s00227-002-0783-z
O. Defeo Æ M. Rueda
Spatial structure, sampling design and abundance estimates in sandy beach macroinfauna: some warnings and new perspectives
Received: 21 September 2001 / Accepted: 10 December 2001 / Published online: 19 February 2002 Springer-Verlag 2002
Abstract We discuss methodological aspects directed to quantify the across-shore population structure and abundance of sandy beach macroinfauna. The reliability of estimates derived from design-based (stratified random sampling) and model-based (geostatistics, kriging) approaches is discussed. Our analysis also addresses potential biases arising from environmentally driven designs that consider a priori fixed strata for sampling macroinfauna, as opposed to species-driven sampling designs, in which the entire range of across-shore distribution is covered. Model-based approaches showed, spatially, highly autocorrelated and persistent structures in two intertidal populations of the Uruguayan coast: the isopod Excirolana armata and the yellow clam Mesodesma mactroides. Both populations presented zonation patterns that ranged from the base of the dunes to upper levels of the subtidal. The Gaussian model consistently explained the spatial distribution of species and population components (clam recruits and adults), with a minor contribution ( £ 5%) of unresolved, smallscale variability. The consistent structure of spatial dependence in annual data strongly suggests an acrossshore-structured process covering close to 35 m. Kriging predictions through cross-validation corroborated the appropriateness of the models fitted through variographic analysis, and the derived abundance estimates were very similar (maximum difference=7%) to those obtained from linear interpolation. Monthly analysis of E. armata data showed marked variations in its zonation and an unstable spatial structure according to the
Communicated by O. Kinne, Oldendorf/Luhe O. Defeo (&) Æ M. Rueda Centro de Investigacio´n y de Estudios Avanzados del IPN, A.P. 73 Cordemex, 97310 Me´rida, Yucata´n, Me´xico E-mail:
[email protected] Fax: +5299-99-812334 O. Defeo UNDECIMAR, Facultad de Ciencias, Igua´ 4225, Montevideo 11400, Uruguay
Gaussian model. The clear spatial structure resulting from species-driven sampling was not observed when data was truncated to simulate an environmentally driven sampling design. In this case, the linear semivariogram indicated a spatial gradient, suggesting that sampling was not performed at the appropriate spatial scale. Further, the cross-validation procedure was not significant, and both density and total abundance were underestimated. We conclude that: (1) geostatistics provides useful additional information about population structure and aids in direct abundance estimation; thus we suggest it as a powerful tool for further applications in the study of sandy beach macroinfauna; and that (2) environmentally driven sampling strategies fail to provide conclusive results about population structure and abundance, and should be avoided in studies of sandy beach populations. This is especially true for microtidal beaches, where unpredictable swash strength precludes a priori stratification through environmental reference points. The need to use adaptive sampling designs and avoid snapshot sampling is also stressed. Methodological implications for the detection of macroecological patterns in sandy beach macroinfauna are also discussed.
Introduction Knowledge about variations in abundance and underlying processes influencing the distribution of organisms in space and time are critical in ecology. Errors in estimating abundance or the spatio-temporal structure of animal populations can induce an erroneous interpretations of and answers to research questions. For invertebrate macroinfauna of sandy beaches, quantitative estimates of abundance are usually straightforward to obtain, because of easy access to beaches and the relatively simple materials used for sampling: no more than a quadrat or a sheet-metal cylinder, a screen and a shovel are needed. Abundance estimates on sandy shores
1216
are usually presented in two ways: (1) as mean densities and (2) as individuals per running meter or strip transect (IST: ind m–1). In the former, the number (or biomass) of individuals per quadrat is averaged by the number of quadrats sampled on the beach, whereas IST estimates the number (or biomass) from each quadrat summed over meters of shoreline by simple linear interpolation, using the following equation (Defeo 1996): n P
qi ISTr ¼ i¼1 wr nr
ð1Þ
where the mean density q (ind m–2 or g m–2) of all n samples pertaining to transect r is multiplied by the corresponding width (w) of the species distribution across the beach. Eq. 1 is used to avoid biased results from using mean density per quadrat (ind m–2) as an abundance estimate, after beach profiles have been impacted during either rough or calm conditions. These variable environmental conditions can result in dramatic contraction/expansion of the across-shore distribution of macroinfauna, implying that density should be considered a biased index of abundance (Brazeiro and Defeo 1996). This emphasizes the need to examine the whole transect, rather than individual sampling quadrats, as the appropriate unit for estimating population abundance (McLachlan 1983; Defeo 1996, 1998; McLachlan et al. 1996). Abundance estimates of macroinfauna need to be expressed as a total quantity, with some measure of the associated variance. The two calculation methods detailed above are based on variance estimates, which assume that the error terms of the samples are stochastically independent of one another (Simard et al. 1992). This approach, which relies on random sampling theory, allows calculations to be done ignoring spatial autocorrelation (small-scale dependence between consecutive sampling units, SUs: Legendre 1993). This condition is not met by sandy beach macroinfauna, which are spatially autocorrelated (Defeo 1993; Gime´nez and Yannicelli 2000). Indeed, most sandy beach populations present strong and persistent (sensu Orensanz and Jamieson 1998) spatial distribution patterns in response to an environment that is spatially and temporally structured by sharp, small-scale gradients. Aggregations persist in time, but, in contrast to sessile species, the relative positions of patches across the beach vary according to the different susceptibilities of each species to variations in environmental conditions (Brazeiro and Defeo 1996, and references therein). In spite of the above, the spatial dimension of population processes and patterns has received little attention in sandy beach ecology, thus overlooking the paradigm of spatial structuring, a critical determinant of how ecosystems function (Legendre 1993). Given the strongly aggregated patterns among sandy beach macroinfauna (Brazeiro 1999a), model-based estimates (e.g. geostatistics) should be relevant when
interpolating abundance between SUs across a beach. Indeed, geostatistics explicitly consider the shape of the spatial autocorrelation in density of consecutive SUs, together with their spatial organization and location on regular grids (surface or volume) or at fixed stations (e.g. single-dimensional transects) (Petitgas 1993; Legendre and Legendre 1998). The fact that mean abundance estimates and their variances are expressed in terms of this structure has made model-based approaches very relevant for crustaceans (Conan 1985; Simard et al. 1992; Maynou et al. 1998), bivalves (see e.g. Thrush et al. 1989; Hewitt et al. 1997) and fishes (Petitgas 1993; Rueda 2001; Rueda and Defeo 2001). Two different sampling designs are usually conducted on sandy beaches: (1) species-driven and (2) environmentally driven. In the former, sampling points are systematically allocated along each of n transects perpendicularly fixed to the shoreline. Transects cover the entire across-shore distribution of the species and at least two or more quadrats in which abundances of zero are sampled in both the upshore and downshore directions before terminating the transect (Defeo 1993, 1998; Go´mez and Defeo 1999; Lercari and Defeo 1999). The second sampling procedure, which is by far the most commonly used, is based on one (Jaramillo et al. 1993) or more (Cardoso and Veloso 1996) shorenormal transects, with SUs randomly allocated at fixed beach levels defined by the position of reference points, such as tidal marks (meso-/macrotidal beaches), swash levels (microtidal beaches), or drift and effluent lines. The validity of this sampling design has been severely questioned (James and Fairweather 1996). Indeed, sandy beach macroinfauna undergoes drastic short-term variations in zonation, associated with the unpredictable intensity of tidal levels generated by onshore winds, storm surges and barometric tides (Brazeiro and Defeo 1996). Considerable changes have also been observed between spring and neap tides (McLachlan 1983). Moreover, sandy beach species tend to occupy different beach levels, according to the season of the year. These short (days) and medium (seasons) term variations in macroinfauna distribution are not always reflecting swash levels or effluent or drift line positions (Brazeiro and Defeo 1996; Gime´nez and Yannicelli 1997). This diminishes the applicability of an environmentally driven sampling design defined by a priori fixed strata, because it does not necessarily cover the entire across-shore species distribution, underestimating the parameter w and, therefore, IST (Eq. 1). The present paper addresses the relative merits of design-based and model-based approaches for estimating abundance in sandy beach populations. This analysis is also discussed in the light of the species- versus environmentally driven sampling designs detailed above. Implications for these methodological approaches in the estimation of population structure and abundance and in the detection of large-scale patterns in sandy beach macroinfauna are also discussed.
1217
Materials and methods Spatial structure and annual abundance estimates To examine the possibility that macroinfauna abundance estimates on sandy beaches may be affected by either the estimation procedure or the sampling design, we used information contained in two databases collected by de Alava and Defeo (1991) and Defeo (1996). In the first case the data was gathered to quantify the intertidal distributional patterns of Excirolana armata (Isopoda: Cirolanidae); sampling was carried out for 1-year (monthly intervals) between May 1988 and April 1989. The second was a longterm study of the yellow clam Mesodesma mactroides, conducted between 1983 and 1990. However, in this paper we used the 12 monthly samplings carried out in 1984 to estimate recruit (43 mm length) densities, averaged for each across-shore level over the entire beach. Both studies were conducted at Barra del Chuy (3340¢S; 5329¢W) in Uruguay, using a systematic sampling design, with 3 (isopods) or 22 (clams) transects perpendicular to the shoreline. No significant differences were ever observed between transects and months, in beach slope, mean grain size, sorting or organic matter of the sediment (two-way ANOVA: P>0.05), so we assume that longshore gradients do not affect across-shore distributions of organisms, as shown elsewhere (Defeo et al. 1986; de Alava and Defeo 1991; Defeo 1993). SUs were spaced every 4 m, from the base of the sand dunes to the sea, until two successive SUs without organisms were recorded. Each SU was sampled with a cylinder (28.2 cm diameter, 40 cm deep), and sediment collected in this manner was sieved through 0.5 mm mesh. All isopods and clams were counted, measured and weighed. The upper (USL) and lower (LSL) limits of the swash zone were determined, respectively, by the maximum and minimum tide advance during sampling time. The across-shore spatial structure of both populations was determined by variographic analysis. Each spatial process consisted of observations measured at a location x defined with respect to its relative position to the upper swash limit (USL=0 m), with negative values indicating a downshore direction. Following the main question addressed in the paper, only one dimension was considered (transects perpendicular to the shoreline), explicitly discarding anisotropic analysis (differential spatial continuity in the across- and longshore directions). We tested for departures from second-order stationarity (sensu Cressie 1991), which require a certain degree of regional homogeneity of abundance and absence of trends (Simard et al. 1992). In this sense, density variations of neither E. armata nor M. mactroides were correlated in the across-shore direction (r between –0.26 and 0.13; P>>0.05). In order to perform variographic analysis, we used a percentage of the maximum lag close to 60% as the default active lag distance. Two reasons induced us to do this: (1) variograms decompose at large lag intervals close to the maximum lag interval (Robertson 2000) and (2) the limited area sampled (average beach width 80 m) and the usually normal distribution of individuals across the dune–sea axis (Defeo et al. 1986, 1997; de Alava and Defeo 1991) produce: (a) few pairs of points available that preclude the obtention of robust estimates and (b) a consistent small variance at large lag intervals, which contradicts the general pattern of constant variance with increasing lags (Fletcher and Sumner 1999). We consistently used a uniform interval of 4 m for variographic analysis, i.e. the minimum distance between two consecutive SUs. A non-directional experimental semivariogram c(h) was estimated by pooling pairs of observations in the downshore direction using Matheron’s (1965) estimator: Nh P
cðhÞ ¼ i¼1
Z^ ðX0 Þ ¼
n X
ki Z ðXi Þ
ð3Þ
i
where ki is the vector of observed densities that minimizes the prediction error (Cressie 1991). Kriging interpolations were evaluated using a jackknife cross-validation. To this end, the observed (O) and estimated (E) densities were plotted and fitted to a linear regression in the form O=a+bE (Cressie 1991). Departures from a 45 line through the origin indicate model inadequacy, and thus the significance of a and b was tested (t-test) under the null hypotheses a=0 and b=1 (Power 1993). Cross-validation was also used for evaluating the appropriateness of the models fitted in the variographic analysis. Abundance of isopods and clams was estimated by: (1) linear interpolation (Eq. 1) and (2) block kriging (Matheron 1965). The latter was used to produce global estimates of average density (ind m–2) and abundance (ind m–1) in discrete areas (blocks) of 1·1 m, over the dune–sea axis. Global abundance estimates were thus obtained through a discrete summation of estimates of interpolated kriging densities (ind m–2) along each entire transect. Standard deviations were obtained following the same reasoning (Simard et al. 1992; Maynou et al. 1998). Zonation dynamics, and environmentally versus species-driven sampling procedures Monthly analyses were performed to test the recurrence of annual patterns. For this purpose we used the database for E. armata, and selected 3 months (May and August 1988 and April 1989) to carry out all calculations detailed for annual data (see ‘‘Results’’ for a detailed explanation of this selection). To compare the spatial structure and abundance estimates between environmentally and species-driven designs, we used E. armata information collected in May 1988, considering: (1) the entire species distribution and (2) only the stations sampled below the USL. This artificially partitioned dataset was used for simulating a priori expected distribution of this typical intertidal species along the Atlantic coast of South America (de Alava and Defeo 1991; Veloso and Cardoso 2001). As Barra del Chuy is a microtidal beach (semidiel tide 0.05), and the correlation coefficient between observed and predicted densities was always highly significant (P
