LJournal of Experimental Marine Biology and Ecology,
216 (1997) 99?128
Spatial structure of bivalves in a sandflat:
Scale and generating processes
,b a a c d*P. Legendre , S.F. Thrush , V.J. Cummings , P.K. Dayton , J. Grant ,
a e f a gJ.E. Hewitt , A.H. Hines , B.H. McArdle , R.D. Pridmore , D.C. Schneider ,
a h aS.J. Turner , R.B. Whitlatch , M.R. Wilkinson
aNational Institute of Water and Atmospheric Research, PO Box 11-115, Hamilton, New Zealand
b
? ? ? ?Departement de sciences biologiques, Universite de Montreal, C.P. 6128, succ. Centre-ville, Montreal,
?Quebec H3C 3J7, Canada
cScripps Institution of Oceanography, UC-SD, La Jolla, CA 92093-0201, USA
dDepartment of Oceanography, Dalhousie University, Halifax, Nova Scotia B3H 4J1, Canada
eSmithsonian Environmental Research Center, P.O. Box 28, Edgewater, MD 21037, USA
fBiostatistics Unit, School of Biological Sciences, University of Auckland, Private Bag, Auckland,
New Zealand
gOcean Sciences Centre, Memorial University, St. John?s, Newfoundland A1C 5S7, Canada
hDepartment of Marine Science, University of Connecticut, Avery Point, Groton, CT 06340-6097, USA
Abstract
A survey was conducted during the summer of 1994 within a fairly homogeneous 12.5 ha area
of sandflat off Wiroa Island, in Manukau Harbour, New Zealand, to identify factors controlling the
spatial distributions of the two dominant bivalves, Macomona liliana Iredale and Austrovenus
stutchburyi (Gray), and to look for evidence of adult?juvenile interactions within and between
species. Most of the large?scale spatial structure detected in the bivalve count variables (two
species, several size classes of each) was explained by the physical and biological variables. The
results of principal component analysis and spatial regression modelling suggest that different
factors are controlling the spatial distributions of adults and juveniles. Larger size classes of both
species displayed significant spatial structure, with physical variables explaining some but not all
of this variation. Smaller organisms were less strongly spatially structured, with virtually all of the
structure explained by physical variables. The physical variables important in the regression
models differed among size classes of a species and between species. Extreme size classes (largest
and smallest) were best explained by the models; physical variables explained from 10% to about
70% of the variation across the study site. Significant residual spatial variability was detected in
the larger bivalves at the scale of the study site. The unexplained variability (20 to 90%) found in
the models is likely to correspond to phenomena operating at smaller scales. Finally, we found no
support for adult?juvenile interactions at the scale of our study site, given our sampling scale,
after controlling for the effects of the available physical variables. This is in contrast to significant
*Corresponding author. Tel.: (514) 343-7591; fax (514) 343-2293; e-mail: Legendre@ere.umontreal.ca
0022-0981/97/$17.00 ? 1997 Elsevier Science B.V. All rights reserved.
PII S0022-0981( 97 )00092-0
100 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
adult?juvenile interactions found in smaller?scale surveys and in field experiments. Our
perception of adult?juvenile interactions thus depends on the scale of study. ? 1997 Elsevier
Science B.V.
Keywords: Adult?juvenile interaction; Autocorrelograms; Austrovenus stutchburyi; Bivalves;
Macomona liliana; Spatial modelling; Spatial structure
1. Introduction
Scale is emerging as one of the critical factors in ecology because our perception of
most ecological variables and processes depends upon the scale at which variables are
measured. A conclusion obtained at one scale may not be valid at another scale without
sufficient knowledge of scaling effects; this is a source of misinterpretation for many
ecological problems (Schneider, 1994). Ecology must deal with scale because organisms
and types of environment are rarely homogeneous. Heterogeneity makes ecological
variables and processes scale-dependent. Environmental forcing, population and com-
munity dynamics, and chance events, are all sources of heterogeneity (Dutilleul and
Legendre, 1993) which contribute to create spatial structures of various kinds, such as
gradients or density-scapes with mountains of high density and valleys of low density
(Schneider, 1987; Legendre and Fortin, 1989; Borcard and Legendre, 1994). The
concept of spatial scale in a sampling design refers to three components: grain, lag and
extent (He et al., 1994; Thrush et al., 1997b). Since field experiments cannot be
conducted at all scales, a good starting point before planning experiments is the
identification of the patterns that can be detected at one or several spatial scales.
In this paper, we assess whether the spatial distributions of infaunal bivalves are
random or spatially structured. If the distributions appear random, for the spatial grain,
lag and extent of the field study, it is unlikely that we will be able to identify
relationships indicating processes important in determining distribution patterns. If,
however, bivalve distributions do exhibit spatial structure, we can formulate hypotheses
about the main determinants of that structure and, by matching the scale of bivalve
spatial variation with those of other variables, provide clues of the scale dependence of
different processes. Correlative studies do not provide conclusive proof of causal
ecological hypotheses, but they may help discard hypotheses for which they provide no
support (assuming the test has enough power). They can also help generate hypotheses
for future experiments to be conducted at the scale(s) suggested by the results.
We formulated a series of hypotheses, from large to small scale, concerning
relationships between bivalve distributions and various factors. Factors may be physical
(large-scale processes, with small to large-scale effects), ranging from variations in tidal
elevation, wind-wave disturbance and tidal current velocity operating at the scale of the
site and beyond, to variation in sediment characteristics around the site (which may also,
ultimately, be controlled by hydrodynamics); or biological (predominantly smaller-scale
processes) such as interactions between species and different post-settlement life-history
stages. The importance of processes may be viewed as corresponding to a gradation in
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 101
effects over decreasing scales, with physical effects predominating at large spatial scales,
while biological effects predominate at smaller scales. Scale-dependent shifts in the
predominance of one process over another (Stommel, 1963; Haury et al., 1978; Amanieu
et al., 1989) have also been referred to as hierarchy theory (Allen and Starr, 1982). Care
must be taken in applying this term to marine systems, however, because they are open
and cannot be described by logical hierarchies. For example, pelagic larvae are not
controlled during their development by the local environmental conditions prevailing at
the locations where they will eventually settle.
In this paper, we examine the spatial distribution of two bivalves, Macomona liliana
Iredale, a deposit- and suspension-feeding tellinid, and Austrovenus stutchburyi (Gray),
a suspension-feeding venerid, within a fairly homogeneous 12.5 ha area of sandflat. In
terms of both density and biomass, Macomona is the most important species on the
Wiroa sandflats. Previous experiments indicate that Macomona plays an important role
in macrofaunal community dynamics (Thrush et al., 1992, 1996) and is an important
food source for eagle rays and waders (Thrush et al., 1994; Cummings et al., 1997;
Hines et al., 1997). We were concerned with determining within rather than between
habitat variation. However, the site incorporated sufficient small-scale variation in
physical features to be representative of the extensive mid-intertidal sandflat habitat of
the region. The extent of the study area and sampling strategy were determined after a
pilot study which examined the spatial scales at which bivalve variability was found in
this habitat (Hewitt et al., 1997). Spatial modelling was used to describe the significant
spatial structures exhibited by Macomona and Austrovenus, assess the consistency of
spatial structures for different size classes, relate patterns to physical factors, and look
for intra- and interspecific relationships between bivalves.
2. Materials and methods
2.1. Study site
A 250 m 3 500 m area (12.5 ha) was selected on the sandflat of Wiroa Island,
Manukau Harbour, New Zealand (378 019 S, 1748 499 E; Fig. 1(a)). A general description
of the area and of its physical characteristics is given in Thrush et al. (1997b) and Bell et
al. (1997).
The area was marked off into 200 grid cells of 25 m 3 25 m each. One sampling
station was selected at random within each cell (with the help of a pseudo-random
number generator) and marked by a peg (Fig. 1(b)). The sampling ??grain?? was a plot of
50 cm 3 50 cm, 15 cm deep, for large bivalves, and three 13 cm diameter 3 15 cm deep
cores for smaller animals. Grain was different again for physical variables, varying from
a point for elevation and for variables derived by modelling, to the cores described
above for sediment composition. The resulting lag (distance between sample centres)
ranged from 5 to 530 m, with a mean distance of 201 m among all pairs of locations,
and a mean distance of 31 m between neighbouring plots.
Sampling was carried out on 22 and 23 January 1994 at the 200 locations described
2
above. On 16 February 1994, new 0.25 m plots were dug out at 31 of the 200 locations
102 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Fig. 1. (a) Position of the study site on the Wiroa Island sandflat. (b) Location of the 200 sampling stations in
the study site, including the 31 stations sampled on 16 February (stars). Coordinates are from an arbitrary zero
mark.
(Fig. 1(b), stars); the rationale through which these 31 locations have been selected is
described in Thrush et al. (1997b). Although tests of significance computed from 31
locations will have lower power than with 200 locations, this information allowed us to
assess the short-term persistence of patterns.
2For each location, three cores totalling 0.04 m were taken and the remaining
2
sediment in the 0.25 m quadrat was excavated to a depth of 15 cm. Sediment collected
in corers was sieved (500 mm mesh) to extract macrofauna, while the remaining
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 103
sediment excavated from the quadrat was sieved (4 mm mesh) to collect large
Macomona liliana and Austrovenus stutchburyi.
2.2. Bivalve count data
For the 22 January data, numbers of Macomona . 15 mm and Austrovenus . 10 mm
are based on the total from three core samples and the remaining sediment excavated
2from the 0.25 m quadrats. Numbers of smaller individuals (Macomona 4?15 mm,
Macomona 2.5?4 mm, Macomona 0.5?2.5 mm, Austrovenus 4?10 mm, Austrovenus
2.5?4 mm, and Austrovenus 0.5?2.5 mm) are based on the three core samples only.
Data were available for 199 of the 200 triplicate core samples only. Data for the missing
sample location were estimated by regression for animals larger than 4 mm, as counts
for quadrats and cores were well correlated for the larger size classes. The 200 core
samples collected in January produced only 17 Austrovenus (4?10 mm), for an
22
estimated mean density of 2.1 animals ? m ; this variable was not analysed.
2On 16 February 1994, the bivalves were counted from 31, 0.25 m quadrats, in the
following categories: Macomona . 15 mm, Macomona 4?15 mm, Austrovenus . 10
mm, and Austrovenus 4?10 mm. This gave us a total of 11 usable bivalve counts: 7 for
22 January and 4 for 16 February. All counts were log-transformed prior to the analyses
(ln(x 1 1)). This transformation was enough to normalise the counts of Macomona . 4
mm and to make all other counts far more symmetrical than the raw data.
2.3. Physical variables
To relate spatial structure in the distribution of Macomona and Austrovenus to various
habitat features, data on a variety of physical variables were collected coincident with
bivalve sampling, with further information on physical processes interpolated from
hydrodynamic models and field measurements (Bell et al., 1997). A major difficulty was
the large number of potentially important controlling variables; in such situations,
choices inevitably have to be made prior to sampling. Our choices were guided by two
weeks of intensive discussions that took place among the authors of this paper and other
invited scientists, during a workshop organised by the National Institute of Water and
Atmospheric Research of New Zealand (NIWA) in Hamilton, NZ, prior to the sampling
campaign itself.
2.3.1. Sediment characteristics
Shell hash (i.e. broken bivalve shell) was measured, at 185 stations, as the dry mass
(g) of broken bivalve shell sieved (500 mm mesh) from the three sediment cores. Shell
hash in our study area was mainly in a layer buried about 3?7 cm below the sediment
surface. Values vary from 5 to 157 g per set of three cores. The shell hash variable was
normalised by the square root transformation. Negative relationships between the shell
hash and adult Macomona could relate to difficulties in moving, and extending siphons
through a shell hash layer. Apparent shell hash effects could also be explained by
covariation with hydrodynamic variables or with elevation.
Data on sediment grain size characteristics (% gravel, % sand and % mud) were
104 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
gathered at 22 locations; only four of these were resampled in February. The study area
is mainly fine sand (98.7%; Thrush et al., 1997b) and previous sampling indicated very
little variation in sediment grain size within the area used for this study (Thrush et al.,
1994). Most if not all the material classified as ??gravel?? for particle size was actually
shell hash.
2.3.2. Elevation
Elevation and position were measured at the 200 grid locations using a geodimeter.
Bed elevation varied by 1.3 m across the 200 sampling stations (from 1.95 to 3.26 m
above chart datum). Elevation is a potentially important variable describing sandflat
topography; it is likely to reflect large-scale zonation patterns.
2.3.3. Hydrodynamics
Hydrodynamic variables are likely to be important determinants of bivalve dis-
tributions influencing larval deposition (Luckenbach, 1984; Butman, 1987; Snelgrove,
1994), the transport of recently settled juveniles (Cummings et al., 1993; Commito et al.,
1995a,b; Roegner et al., 1995), food supply (Emerson and Grant, 1992), and feeding
behaviour (Ertman and Jumars, 1988; Monismith et al., 1990; O?Riordan et al., 1993).
As an important determinant of sediment grain size characteristics, hydrodynamic forces
are likely to covary with sediment variables (Snelgrove and Butman, 1994). The
variables that were available to the present study were derived from numerical model
simulations for tidal currents and wind-waves, described by Bell et al. (1997). Two
predictions can be made about physical?biological interactions: (1) more juveniles
should be found where bed shear stress and wave action are lower; (2) regions of higher
bed shear stress and wave action may be preferred by adults, because physical energy
may maximise the supply of food.
Waves and currents both generate shear stress at the sediment surface (Grant and
Madsen, 1979), but for convenience we treat the two processes separately. The drag
force exerted upon water moving over the sea bed demands that the moving fluid impart
some of its momentum to the seabed. At the sediment interface, the transfer of
momentum (i.e. shear stress) is a maximum, and in turbulent flows it is proportional to
2 2the square of the time-averaged fluid velocity (U ), with the stress equal to rC Ux
(where C is the drag coefficient and r is seawater density). For tidal currents, bed shearx
22
stresses (N ? m ) under peak ebb- and flood-tide velocities, during a mean tide,
correspond to shearing forces applied per unit area; they were computed from the
depth-averaged tidal hydrodynamic model for a mean tide. A further variable was
computed in the form of rate of energy dissipation per unit area at the bed, or the power
3 23
expended per unit area (proportional to U ); it is measured in kg ? s . Both variables
have been multiplied by 1000 for convenience. While the distributions of the peak ebb
tide variables were fairly normal-looking, this is not the case for the peak flood tide
variables (shear stress and rate of energy dissipation) that were skewed positively; a log
transformation solved this problem (ln(x)).
Wind-waves also disturb the sandflat and generate sediment transport. For waves, the
rate of energy dissipation has been integrated over time during a tidal cycle when the
22bed was inundated, to give work done per unit area of bed (kg ? s ). The drag
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 105
coefficient is typically an order of magnitude higher under waves than for tidal currents.
Values were calculated based on model simulations of peak wave orbital velocities at the
bed condensed to a single scenario of a 0.3 m wave height (exceeded only 20% of the
time at this site for onshore winds $ 17 knots) and a typical mean period of 2.7 s. This
was done for the two most common wave directions, SW and WSW (Bell et al., 1997).
The percentage of time each station of the grid is covered by spring and neap tides was
also computed, as well as the percentage of time the waves stir the plot during spring
and neap tides (R.A. Walters, pers. comm.). All these variables are strongly and
negatively correlated with elevation, as expected.
Only some of the physical variables were included in the regression analyses. Among
the 22 physical variables mentioned above, many are likely to be highly correlated. For
each group of physical variables, those that were more strongly correlated with bivalve
counts were selected for inclusion in statistical modelling. Preventing different but
correlated variables from becoming significant in different models makes the comparison
of models easier. The pre-selection procedure presents the same problems as forward
selection of explanatory variables prior to modelling. Hopefully, this risk will be
counterbalanced by a gain in clarity of the models. The alternative would have been to
resort to ridge regression to deal with high collinearities; this methods presents problems
of its own (Legendre and Legendre, 1997).
1. Six highly collinear ??water cover?? variables were available. The ??percent of time the
plot is covered by more than 20 cm water during spring tide?? had the largest sum of
correlations with the other five, as well as the largest correlations (when significant)
with all the bivalve count variables. This variable is biologically reasonable,
reflecting the time when food particles are put into motion by wave activity. This
variable only was used in the statistical model.
2. Three highly collinear ??wave stirring?? variables were available. Among them, the
??percent of time large waves stir the plot during spring tide?? had the largest sum of
correlations with all the others. It is also strongly collinear with the water cover
variable retained in the previous paragraph (r 5 0.953).
3. Peak ebb shear stress and rate of energy dissipation were very highly correlated
(r 5 0.9995), as were peak flood shear stress and rate of energy dissipation (r 5
0.9999). Thus, only the shear stress variables, that relate directly to the potential for
sediment transport in the absence of waves, were retained in the spatial modelling, in
order to reduce collinearity.
4. The work variables (SW and WSW winds) were retained in the spatial modelling.
Correlations with the bivalve counts were high.
The preliminary analysis, which started with 22 physical variables, resulted in eight
variables that were spatially structured (they all had highly significant trend surface
equations: Section 2.4) and significantly related to the bivalve count data. Elevation to
the powers 2 and 3 were also included in the modelling effort in order to model
nonlinear relationships to bivalve counts. The correlations among the variables to be
used in modelling are summarised in Table 1. The large correlations (r . 0.6),
underscored in that table, indicate that there is still a large amount of collinearity among
106
P
.Legendre
et
al
.
/
J
.Exp
.M
ar
.Biol
.Ecol
.216
(1997)99
?128
Table 1
Pearson correlations among the physical variables retained for spatial modelling of the bivalve count data
2 3Shell hash Elevation Elevation Elevation Ebb stress Flood stress SW work WSWwork .20 cm water Wave stirring
aShell hash 1.0000 0.3505 0.3404 0.3286 0.2784 20.4935 0.3842 0.2912 20.3417 20.4050
Elevation 0.3505 1.0000 0.9980 0.9926 0.1266 20.6500 0.9002 0.8985 20.9992 20.9556]] ]] ]] ]] ]] ]] ]]2Elevation 0.3404 0.9980 1.0000 0.9983 0.1381 20.6588 0.8966 0.9211 20.9997 20.9547]] ]] ]] ]] ]] ]] ]]3Elevation 0.3286 0.9926 0.9983 1.0000 0.1470 20.6649 0.8893 0.9392 20.9966 20.9495]] ]] ]] ]] ]] ]] ]]
Ebb stress 0.2784 0.1266 0.1381 0.1470 1.0000 20.4092 0.0365 0.1975 20.1309 20.2235
bFlood stress 20.4935 20.6500 20.6588 20.6649 20.4092 1.0000 20.4864 20.7243 0.6566 0.6084]] ]] ]] ]] ]] ]]
SW work 0.3842 0.9002 0.8966 0.8893 0.0365 20.4864 1.0000 0.7795 20.8966 20.8782]] ]] ]] ]] ]] ]]
WSW work 0.2912 0.8985 0.9211 0.9392 0.1975 20.7243 0.7795 1.0000 20.9139 20.8507]] ]] ]] ]] ]] ]] ]]c.20cm water 20.3417 20.9992 20.9997 20.9966 20.1309 0.6566 20.8966 20.9139 1.0000 0.9531]] ]] ]] ]] ]] ]] ]]dWave stirring 20.4050 20.9556 20.9547 20.9495 20.2235 0.6084 20.8782 20.8507 0.9531 1.0000]] ]] ]] ]] ]] ]]
n5200, except in comparisons involving shell hash (n5185). Correlations larger than 0.6 are underscored.
a Shell hash: square root transformation. n5185.
b Flood shear stress: natural logarithm transformation.
C % time the plot is covered by more than 20 cm water during spring tide.
d % time large waves stir the plot during spring tide.
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 107
these variables; only shell hash and ebb shear stress are fairly linearly independent of the
other physical variables. Collinearity was reduced by backward elimination of non-
significant variables during regression modelling.
2.3.4. Chemistry
To indicate the availability of food to deposit feeders, data were obtained at 19
stations for percentage of organic carbon and nitrogen in the sediment, using a
Perkin?Elmer elemental analyser. Correlations between these variables and log-trans-
formed bivalve counts were in general low. These variables, which are not available for
the 200 stations, were not used in the modelling. Surficial sediment chlorophyll
concentration was also determined following the methods described in Thrush et al.
(1994), but after analysing a random subset, no variation had been detected and the
variable was abandoned.
2.4. Statistical methods
A variety of methods have been proposed to investigate spatial autocorrelation,
including spectral analysis, spatial autocorrelograms, variograms and other forms of
variance?distance curves; these techniques are presented in various textbooks and
review papers about spatial statistics and numerical ecology (Sokal and Oden, 1978;
Cliff and Ord, 1981; Ripley, 1981; Upton and Fingleton, 1985; Burrough, 1987;
Legendre and Fortin, 1989; Isaaks and Srivastava, 1989; Haining, 1990; Legendre and
Legendre, 1997). Frequency-based techniques (e.g. spectral analysis) cannot be used for
irregularly spaced observations; distance-based techniques (e.g. correlogram analysis)
are more appropriate here.
A combination of spatial autocorrelograms, trend surface analysis, and mapping was
used to describe spatial structures. Spatial autocorrelograms were computed using
Moran?s I spatial autocorrelation coefficient; following Oden (1984), an overall test of
the significance of each spatial autocorrelogram was performed by checking that the
most significant spatial autocorrelation coefficient found in a correlogram was significant
at a Bonferroni-corrected significance level a9 5 a /k where k is the number of
autocorrelation coefficients in the correlogram.
Trend surface analysis (Student, 1914) is a regression of a dependent variable y on a
polynomial function of the geographic coordinates X and Y of the sampling stations
where the variable has been observed or measured. The X and Y coordinates were
easting and northing respectively, in km, measured from an arbitrary ??zero?? surveyor?s
mark. These coordinates were centred on their respective means in order to reduce the
linear dependency (collinearity) between the first and second-degree terms of the spatial
polynomial of geographic coordinates; the mean point corresponds roughly to the centre
of the study site. The amount of variation explained by a trend surface equation is not
changed by a translation of the spatial coordinates across the map. The trend surfaces
were computed by ordinary least-squares fitting of a polynomial equation of these
centred X?s and Y?s. The following third-degree polynomial equation was used in the
present study, as suggested by Legendre (1990):
2 2 3 2 2 3y? 5 b 1 b X 1 b Y 1 b X 1 b XY 1 b Y 1 b X 1 b X Y 1 b XY 1 b Y0 1 2 3 4 5 6 7 8 9
108 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
where the b?s are the regression coefficients to be estimated by regression. Non-
significant terms (called monomials) of the spatial polynomial were removed using a
backward elimination procedure.
The rationale behind this polynomial is the following. Let us assume some general
shape for the biological phenomena to be described; for instance, a phenomenon may
start from some mean value of the measured variable, increase in intensity to a
maximum, then go down to a minimum, and come back to the mean value. The amount
of space required for the phenomenon to complete a full cycle?whatever the shape it
may take?is called its scale. Commonly used models for such shapes are sine functions;
models for these functions are easy to generate, and several methods exist in the time
series literature for fitting them to actual data series. In the present study, we will not
restrict allowable phenomena to trigonometric functions; we will try to model instead
any phenomenon that has the general behaviour described above (mean, maximum,
minimum, and back to the mean) using polynomial functions; these functions are more
flexible than sines or cosines, in that they do not require symmetry or strict periodicity.
The degree of the polynomial which is appropriate to model an anticipated phenomenon
is predictable. For instance, if a variable has spatial variation at the scale of the study
site (say, in the X direction5easting), it should be correctly modelled by a polynomial of
degree 3, which has two extreme values, a minimum and a maximum. If the scale of the
phenomenon is larger than the study site, a polynomial of degree less than 3 should be
sufficient; degree 2 if only one maximum, or only one minimum, is observed in the
sampling window; and degree 1 if the study site is limited to the increasing, or
decreasing, portion of the phenomenon. Conversely, if the scale of the phenomenon
controlling the variable is smaller than the study site, more than two extreme values
(minima and maxima) should be found in the study site, and a polynomial of order larger
than 3 would be required to model it correctly. So, using a polynomial of degree 3 acts
as a filter, because it is a way of looking for phenomena that are of the same scale, or
larger, than the study site. The same reasoning applies to the X (5easting) and Y
(5northing) directions if we use of a polynomial combining the X and Y geographic
coordinates.
Correlations among bivalves in different size classes were investigated using principal
component analysis. Eigenvalues and eigenvectors were computed from the correlation
matrix and the eigenvectors were scaled to the square root of their respective
eigenvalues. With this scaling, plots correctly represent the projection angles among
variables. Arc cosine transformations of the correlations give the angles between
variables in multidimensional space (Legendre and Legendre, 1983). Visually, if lines
are drawn from the origin to the points representing different variables, then the angle
between two lines represents the correlation between the corresponding variables, a
small angle meaning a high correlation.
Spatial modelling was performed using a method derived from that proposed by
Borcard et al. (1992) and Borcard and Legendre (1994) to partition the variance of a
dependent variable (or set of dependent variables) among environmental and spatial
components. Our specific purpose here was to take the spatially-structured variation as a
measure of the dependent variables? variation worth explaining, and to measure by how
much this fraction had decreased after incorporating different groups of explanatory
variables into the models.
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 109
2.5. Spatial modelling
Each bivalve count variable was modelled through a process involving five steps. The
physical variables were used first to explain the spatial variation of the bivalve counts,
followed by the biological variables. For each dependent bivalve count variable, only
those size classes that were larger than or equal to that of the modelled variable were
used as explanatory variables, following the hypothesis that small animals can only be
directly affected by other bivalves of the same or larger sizes (e.g. feeding adults may
disturb juveniles but not the converse). Previous studies on this sandflat provided no
indication for indirect effects of aggregations of juveniles on adults. Biological
interaction variables may also explain part of the small-scale, non-spatially structured
variation in the data. When biological interactions were found, we checked whether the
variation thus explained was spatial or not.
21. Step 1: Calculate the coefficient of multiple determination (R ) of the spatially
structured variation in the dependent variable.
2. Step 2: Do the multiple regression modelling, selecting from the physical variables
described in the previous section.
3. Step 3: Check whether the remaining variation is spatially structured, by adding the
nine terms of the spatial polynomial to the multiple regression. The statistic of
2interest is the increase in explained variation, DR . When a significant spatial
2polynomial cannot be found, the DR value is given for all nine terms of the spatial
polynomial.
4. Step 4: Add the biological interaction variables to the physical equation, if significant.
2Two statistics are of interest: the total fraction of variation, R , explained by the
combined action of the physical and biological variables, and the increase in
2
explained variation due to the biological variables alone DR (with indication of the
significance of that increase).
5. Step 5: If the biology provided significant explanatory variables, does there remain
some spatially-structured variation in the residuals? If so, it could correspond to
spatial variation not explained by the variables in the model, or to large-scale spatial
variation appearing in the residual data after removing local biological effects. The
2
statistic provided, DR , is the same as in step 3.
6. Summary statistic: total explained variation at the outset of the modelling process
2(R ), for significant variables, excluding the variation explained by non-significant
spatial polynomials.
3. Results
3.1. Do the distributions of bivalves display significant spatial structures?
3.1.1. Spatial autocorrelograms of bivalve counts
The spatial all-directional autocorrelograms (log-transformed data) for bivalve counts
are presented in Fig. 2. All correlograms for 22 January, based upon a large number of
observations, were globally significant. In contrast, the correlograms for 16 February
110 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Fig. 2. All-directional spatial autocorrelograms for log-transformed bivalve counts. Graphs (a?d) are based
upon 199 or 200 sites sampled on 22 January, while graphs (e?f) are based upon 31 sites sampled on 16
February. Significant values of Moran?s I spatial autocorrelation coefficient ( p#0.05) are represented by black
symbols.
were generally not significant, except the one for Macomona.15 mm. Lack of
significance can be attributed to low power, due to the calculations being based upon
only 31 sampling stations. Recomputing these correlograms with 10 distance classes
instead of 20 did not lead to more significant results.
Interpretation of the spatial structures represented by the significant correlograms was
based on the simulations presented in Legendre and Fortin (1989). Fig. 2(a) and 2(b)
correspond to large aggregation structures (shaped as bumps, troughs, or waves) that are
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 111
Table 2
Trend surface models for the bivalve count variables (ln(x11))
2d Model for Macomona.15 mm, 22 January, counted in 0.25 m
2
?n5200, o58003, y540.02, S.E.50.791, dens.5160.1 R 50.3171, p,0.0001
3 2y?53.69421.927 X11.722 Y27.590 XY136.230 X 247.336 X Y
2d Model for Macomona 4?15 mm, 22 January, counted in 0.04 m
2
?n5200, o5209, y51.05, S.E.50.085, dens.526.1 R 50.1192, p,0.0001
2y?50.68410.671 X11.319 Y26.161 X
2d Model for Macomona 2.5?4 mm, 22 January, counted in 0.04 m
2
?n5199, o5139, y50.70, S.E.50.080, dens.517.5 R 50.1094, p,0.0001
2y?50.42910.889 X24.288 X 110.135 XY
2d Model for Macomona 0.5?2.5 mm, 22 January, counted in 0.04 m
2
?n5199, o5551, y52.77, S.E.50.207, dens.569.2 R 50.1458, p,0.0001
2y?51.14710.924 X11.562 Y29.212 X 119.395 XY
2d Model for Macomona.15 mm, 16 February, counted in 0.25 m
2
?n531, o5964, y531.10, S.E.52.162, dens.5124.4 R 50.3936, p50.0033
2y?53.42420.939 X14.012 Y2133.062 X Y
2d Model for Macomona 4?15 mm, 16 February, counted in 0.25 m
2
?n531, o5202, y56.52, S.E.50.753, dens.526.1 R 50.1503, p50.0312
2y?51.467115.412 X
2d Model for Austrovenus.10 mm, 22 January, counted in 0.25 m
2
?n5200, o5374, y51.87, S.E.50.291, dens.57.5 R 50.3483, p,0.0001
2 2 3 2y?50.28811.630 X118.088 X 231.653 XY123.727 Y 249.721 X 187.807 X Y
2d Model for Austrovenus 2.5?4 mm, 22 January, counted in 0.04 m
2
?n5199, o5185, y50.93, S.E.50.084, dens.523.2 R 50.1617, p,0.0001
2y?50.49521.210 X13.127 Y289.791 X Y
2d Model for Austrovenus 0.5?2.5 mm, 22 January, counted in 0.04 m
2
?n5199, o5192, y50.96, S.E.50.095, dens.524.1 R 50.1694, p,0.0001
2 2y?50.48221.718 X19.415 XY267.507 X Y1131.171 XY
2d Model for Austrovenus.10 mm, 16 February, counted in 0.25 m
2
?n531, o572, y52.32, S.E.51.131, dens.59.3 R 50.5508, p,0.0001
2y?50.119127.920 X 240.770 XY
2d Model for Austrovenus 4?10 mm, 16 February, counted in 0.25 m
2
?n531, o556, y51.81, S.E.50.351, dens.57.2 R 50.2917, p50.0080
3y?50.75619.493 Y2702.473 Y
X and Y are geographic coordinates in km, centred on their respective means. The fitted value of the dependent
2
?variable in each model is designated by y. R is the coefficient of multiple determination of the model, and p
the associated probability. All terms reported in the trend surface equations are significant ( p#0.05); some
regression coefficients are large because the squared and cubic terms of the spatial polynomial are very small
numbers. Basic statistics are also provided for the untransformed counts: n is the number of observations; o
?designates the sum of bivalves of the given species and size, counted in all samples; y is the mean, and S.E. is
22the standard error of the mean; dens. is the mean estimated density (animal?m ).
the same size or larger than the study site; for bumpy or wavy structures, the distance
between successive peaks or troughs is twice the distance where the minimum
autocorrelation value occurs in the correlogram. Fig. 2(c) and 2(d) correspond to spatial
gradients running over the site, that are perhaps part of bumpy or wavy structures that
occur on spatial scales larger than our study area. Correlograms cannot be of further help
in determining the shapes of these large-scale structures because they exceed the size of
the study site. So, their shapes will now be investigated using trend surface analysis and
mapping; Figs. 4?7 will help interpret the correlograms.
112 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
2Fig. 3. Coefficients of determination (R ) of the spatial trend-surface models, where a significant model has
been identified. Details in Table 1.
3.1.2. Trend surface analyses
Highly significant trend surface equations were found for all bivalve variables (Table
2). This corroborates the interpretations of the autocorrelograms; the spatial distributions
of these organisms are not random, but highly organised at the scale of the 12.5 ha study
site. The trend surface models for the smaller animals have much smaller coefficients of
determination (10?20%) than for larger animals (30?55%). The best models, that is, the
2
models with the highest coefficients of determination (R ), are for the Macomona.15
mm and Austrovenus.10 mm (Table 2, Fig. 3). Also, the coefficients of determination
are consistently higher for Austrovenus than for Macomona, despite the fact that
Macomona were usually far more numerous than Austrovenus. That the coefficients of
determination are consistently higher for 16 February than for 22 January is due to a
large extent to a higher ratio of number of parameters to sample size.
3.1.3. Maps
Maps that illustrate the trend surface equations are presented for the largest and
smallest size classes (Fig. 4 Fig. 5 Fig. 6 Fig. 7); the field counts are also presented in
each case for comparison. The two species displayed very different spatial patterns for
their largest size classes (the correlation of these two variables across 200 locations is
low, r50.0481), although they presented very similar correlograms (Fig. 2 (a?b));
Legendre and Fortin (1989) had already shown that different spatial structures may lead
to the same type of correlogram. What the two species had in common is a wavy
structure whose main axis of variation is NW to SE. In the same way, while the trend
surfaces for the smallest Macomona and Austrovenus are very different (the correlation
of these two variables is r50.2039), the all-directional correlograms seem to have
picked up mostly the main spatial gradient in each of these surfaces (Fig. 2 (c?d)).
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 113
Fig. 4. Macomona.15 mm from 200 sites, 22 January. (a) Actual counts at sampling stations in the 200
regular grid cells; in the field, the stations were not equispaced. (b) Map of the trend surface equation
explaining 32% of the spatial variability of the data. The values estimated from the trend-surface equation
(log-transformed data) have been back-transformed to raw counts before plotting. The sampling grid is viewed
from the south.
3.2. Are the patterns stable through time?
Comparisons between 22 January and 16 February are possible only for large animals
(Macomona.15 mm, Austrovenus.10 mm), because only these variables provide
whole-plot counts for both dates, even if only at 31 stations. Paired t-tests performed on
the normalised (log-transformed: ln(x11)) data show that there is a slight but highly
significant difference in means of the log-transformed data for large Macomona
(abundances decreased: mean for 22 January in the original count scale540 animals per
quadrat, mean for 16 February531; for log-transformed data: t55.27, d.f.530; the
difference would remain significant after correction for autocorrelation), but not for large
Austrovenus (mean for 22 January51.9 animals per quadrat, mean for 16 February52.3;
for log-transformed data: t50.81, d.f.530).
114 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Fig. 5. Macomona 0.5?2.5 mm from 199 sites, 22 January. Presentation as in Fig. 3. Dot: missing datum. The
trend surface equation explains 15% of the spatial variability of the data.
To assess the consistency of spatial patterns, for the 31 stations that were sampled on
two occasions, a spatial pattern of differences was determined (Legendre and McArdle,
1997). For both large Macomona and large Austrovenus, no significant trend surface
could be identified among the 31 difference values. So there is no indication of changes
in shape over the three-week interval. This exercise could not be done for the next size
class as data were not available for the same sampling grain (see Bivalve count data,
Section 2.2).
Correlations between dates were also computed to measure the similarity, or stability
of values between dates; differences in mean would not be found by this analysis
because correlations are computed on differences from the respective means. Correla-
tions were high for large Macomona, between whole-plot counts on 22 January and 16
February (r50.7040); the same was true for large Austrovenus (r50.7050). So, again,
adult counts seem fairly stable across three weeks. Such was not the case for smaller
individuals, however; the correlation for Macomona 4?15 mm between core counts on
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 115
Fig. 6. Austrovenus.10 mm from 200 sites, 22 January. Presentation as in Fig. 3. The trend surface equation
explains 35% of the spatial variability of the data.
22 January and whole-plot counts on 16 February is negative (r520.3055). This
indicates that a change occurred.
3.3. Are the patterns the same across species and size classes?
The principal component analysis (Fig. 8(a)) clearly reveals that large Macomona
(.15 mm) had similar distributions in the two dates; the same is true for large
Austrovenus (.10 mm). However, the large Macomona behave very differently from
the large Austrovenus. The same picture was obtained using log-transformed data instead
of raw counts.
Drawing arrows from the smaller to the larger size classes of each species, for each
sampling date, suggests an interesting relationship (Fig. 8(b?c)). While intermediate-
sized Macomona (4?15 mm) on 22 January are unrelated to intermediate-sized
Macomona on 16 February, the large animals are quite highly correlated. Because there
were so few Austrovenus 4?10 mm in the core samples of 22 January, we don?t know
116 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Fig. 7. Austrovenus 0.5?2.5 mm from 199 sites, 22 January. Presentation as in Fig. 3. Dot: missing datum. The
trend surface equation explains 17% of the spatial variability of the data.
whether animals pertaining to this size class would be found in the same or different
portions of the graph on 22 January and 16 February, although we know that their mean
22density over the study site changed from 2 to 7 animal?m during that time interval. We
know, however, that large Austrovenus (.10 mm) are clearly found together in the plot
on the two dates, in a location opposite to large Macomona. In contrast, small animals
(,4 mm) of both species are found in the same region of the plane of the first two
principal components (lower left in Fig. 8(a)); their correlation is 0.2039.
3.4. Spatial modelling of bivalve counts
A summary of the five-step spatial modelling procedure, for each of the 11 bivalve
count variables, is presented in Table 3, with detailed examples in Appendix A.
Interpretation will focus on the signs of the significant regression coefficients, summa-
rised in Table 4. Besides the variables found in Table 4, variables from the spatial
polynomial turned out to be significant in three models (modelling step 3 in Table 3;
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 117
Fig. 8. (a) Bivalve count variables represented in space of the first two eigenvectors, accounting together for
40% of the variation of the correlation matrix (1: 20.4%; 2: 19.6%). Macomona is represented by circles,
Austrovenus by squares. (b) Interpretation of the same graph for Macomona. Arrows indicate the size sequence
for each sampling date. Size classes are represented by numbers. (c) Same for Austrovenus.
only the signs of the significant partial regression coefficients are shown): 2X, 2XY and
3 21X for Macomona.15 mm, 22 January; 1X for Macomona 4?15 mm, 16 February;
3 22X for Austrovenus.10 mm, 22 January; 1X and 2XY for Austrovenus.10 mm,
16 February.
At the scale of the 250 m3500 m study site, the physical variables had significant
contributions to the explanation of count variation in all Macomona and Austrovenus
size classes, except Macomona (4?15 mm) sampled on 16 February (Table 3). The
biological variables explained residual variation in counts in the smallest Macomona and
Austrovenus.
Residual spatial variability, operating at the scale of the study site or larger, was
detected in the larger animals. In the largest Macomona and Austrovenus, for example,
3
spatial terms in X were significant in modelling step 3, indicating the presence of
residual spatial variation at the scale of 500 m. Combinations of physical and biological
variables explained most of the large-scale spatial structure in bivalve counts (scale of
the study site); indeed, there is little or no significant spatial variation left in modelling
2
steps 3 and 5 of Table 3. The unexplained variability (1 minus the R value given in last
column of Table 3) can be attributed either to spatially-structured variability at the scale
of the study site that cannot be expressed as a linear combination of the terms of the
118 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Table 3
Spatial modelling of the bivalve count variables (dependent variables)
Step 1 Step 2 Step 3 Step 4 Step 5
Spatial Physics 1spatial Physics1biology 1spatial Total
2 2 2 2 2 2 2Dependent variable n R R DR R DR DR R
Macomona
.15 mm, 22 Jan. 200 0.3171 0.2583* 0.0726* ? ? ? 0.3309
4?15 mm, 22 Jan. 200 0.1192 0.1060* 0.0268ns ? ? ? 0.1060
2.5?4 mm, 22 Jan. 199 0.1094 0.1653* 0.0345ns ? ? ? 0.1653
0.5?2.5 mm, 22 Jan. 199 0.1458 0.2234* 0.0233ns 0.3785* 0.1551* 0.0320ns 0.3785
.15 mm, 16 Feb. 31 0.3936 0.4880* 0.0926ns ? ? ? 0.4880
4?15 mm, 16 Feb. 31 0.1503 ? 0.1503* ? ? ? 0.1503
Austrovenus
.10 mm, 22 Jan. 200 0.3483 0.3421* 0.0186* ? ? ? 0.3606
2.5?4 mm, 22 Jan. 199 0.1617 0.2126* 0.0304ns ? ? ? 0.2126
0.5?2.5 mm, 22 Jan. 199 0.1694 0.1894* 0.0316ns 0.2810* 0.0917* 0.0378ns 0.2810
.10 mm, 16 Feb. 31 0.5508 0.6872* 0.1036* ? ? ? 0.7907
4?10 mm, 16 Feb. 31 0.2917 0.3437* 0.2383ns ? ? ? 0.3437
2Step 1 reports the R coefficients of the spatial models in Table 2. Step 2: model using the physical variables
2 2
only: R . Step 3: adding the spatial to the physical variables: DR is reported; when a significant spatial
2polynomial cannot be found, the DR is given for all nine terms of the spatial polynomial. Step 4: the physical
2 2 2
and biological variables: R , DR . Step 5: adding the spatial to the physical and biological variables: DR .
2 2Total R : R value reached using all the significant variables in the models. The significance level is 0.05.
2 2 2R 5coefficient of determination; DR 5increase in R from the model without the stated ??1?? variables to the
model that includes them. *5significant at the 0.05 level; ns5spatial structure not significant; ?5no
significant term was found.
cubic trend-surface equation, or to smaller-scale phenomena, since the spatial polyno-
mial of the geographic coordinates X and Y was limited by design to power 3, as
explained in Section 2.4; smaller scales have been investigated by Hewitt et al. (1997).
4. Discussion
We have shown (Fig. 2; Figs. 4?7) that the bivalve count variables are spatially
structured at the scale of our study site. We expected physical processes to be also
spatially structured at that scale, because of the large size of the study site. So, the action
of physical variables on bivalve counts, if any, should be detectable at that scale, and
could contribute to explain the spatial structure of the bivalve counts, detected by the
trend-surface polynomial equations (Table 2).
Our pilot study (Hewitt et al., 1997: cross-correlograms) indicated that interactions
between species and size classes of the two dominant bivalves, when they could be
detected, were found at scales no larger than 0 to 5 m. Negative adult?juvenile
interactions amongst Macomona were indicated by the increased flux of juveniles in
areas of high adult density (Turner et al., 1997) and were demonstrated in experimental
plots throughout the study site (Thrush et al., 1997a). Thus interactions appear to
P
.Legendre
et
al
.
/
J
.Exp
.M
ar
.Biol
.Ecol
.216
(1997)99
?128
119
Table 4
Summary of the significant contributions of the physical and biological variables to the various regression models
2 3Shell Elevation Elevation Elevation Ebb Flood SW WSW .20 cm Wave Macomona 2.5?4 mm Maconoma Austrovenus Austrovenus
hash stress stress work work water stirring 0.5?2.5 mm 2.5?4 mm 0.5?2.5 mm
Macomona
.15 mm, 22 Jan. 1 2 1 1
4?15 mm, 22 Jan. 1 1
2.5?4 mm, 22 Jan. 2 1 1 1
0.5?2.5 mm, 22 Jan. 2 1 1 1 1 1 1 1
.15 mm, 16 Feb. 1 1 1
4?15 mm, 16 Feb.
Austrovenus
.10 mm, 22 Jan. 2 1 2 1 2
2.5?4 mm, 22 Jan. 2 1 2 1 1 1
0.5?2.5 mm, 22 Jan. 1 2 1 1 2 2 1
.10 mm, 16 Feb. 2 1 2 2
4?10 mm, 16 Feb. 2 2
Signs are the signs of the partial regression coefficients.
120 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
function well below the scale of the present study, where the mean distance among
neighbouring samples is about 30 m. As a consequence, we do not expect biological
interactions between the life stages of the two dominant bivalve species to explain any
of the observed large-scale variation; this point is further discussed in Section 4.2. On
the other hand, if biological variables are incorporated into models as surrogates for
unmeasured physical processes (e.g. localised phenomena not effectively modelled or
measured by Bell et al., 1997), they should not create any new large-scale spatial
structure in the residuals of the models. This point has been examined (modelling step
5).
4.1. Explaining the spatial structure of bivalve counts
Spatial variation (i.e. the variation explained by the spatial polynomial) is used in the
present study simply as an indication of a significant spatial pattern in the dependent
variable under study (bivalve count), at the scale of our study site. Such an indication
does not contain any interpretation per se. It simply justifies our search for ecologically
meaningful hypotheses capable of explaining away the spatially-structured variation of
the dependent variable (Borcard and Legendre, 1994). The spatial modelling is
considered fully successful when there is no significant spatial variation left to be
explained in the data. Admittedly, the technique is limited in that the spatial polynomial
used in the present study can only capture the large-scale structures of the dependent
variables. Other techniques should be used to model small-scale autocorrelation in the
data (Legendre and Borcard, 1994).
Our analyses have shown that there are large-scale spatial structures in all bivalve
count variables investigated in this study (Table 2; Table 3, step 1), and that most of
these large-scale structures disappear when physical and biological variables are
included in models. This result also indicates that our sampling design?25 m resolution
within a 250 m by 500 m study site?was appropriate because effects of physical
variables were detected on bivalves sampled at that scale.
Significant large-scale spatial variation remained unexplained by our models for large
bivalves, even though we had included variables for all of the physical factors
commonly evaluated by benthic ecologists. Patterns in the maps of the model residuals
(Fig. 9 Fig. 10) may indicate possible explanations for the remaining variation (Borcard
and Legendre, 1994). In the present case, the trend surface equations of the residuals
3indicate that some process, structured as X , seems to be at work in both maps, but with
effects bearing opposite signs. The nature of that process remains unknown.
The importance of ecological processes (physical or biological) varied for different-
sized animals, across a study site chosen to represent basically a flat and homogeneous
environment in terms of its major physical characteristics. (a) Larger organisms
displayed a significant spatial structure, with physical variables explaining some but not
all of this spatial variation. (b) Smaller organisms were less strongly spatially structured;
virtually all of their large-scale spatial variation could be explained by physical
variables. (c) Regression models for larger organisms include physical variables only,
whereas models for the smallest animals include physical variables as well as a few
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 121
3Fig. 9. Map of the spatial component model (20.920 X27.914 XY123.854 X ) of Macomona.15 mm, 22
2January. It corresponds to DR in step 3 of Table 3. Compare to Fig. 4.
biological variables (other juveniles). The physical variables that are included in the
regression models differ among size classes of a species and between species.
4.2. Biological interactions
Interspecific interactions did not explain any of the spatial variation of the larger
bivalves. The smallest size classes only (0.5?2.5 mm) responded significantly to inter-
and intraspecific interactions amongst juveniles (Table 4). The significant interactions
were positive, except between the smallest Austrovenus (0.5?2.5 mm) and slightly larger
Macomona (2.5?4 mm). No avoidance of adults on the part of juveniles has been
detected in these models, after the physical variables had been taken into account.
3Fig. 10. Map of the spatial component model (252.745 X ) of Austrovenus.10 mm, 22 January. It
2
corresponds to DR in step 3 of Table 3. The slight increase along the 02250 m axis is because the sampling
grid is not aligned with the easting (X) and northing (Y) geographic coordinates (Fig. 1). Compare to Fig. 6.
122 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Mechanisms for juvenile?juvenile interactions have been documented in some labora-
tory studies (Ahn et al., 1993), but these relationships may not mean that biological
interactions actually took place; they may simply indicate that there were other
unmeasured physical variables operating at this small scale, that caused the juveniles to
deposit and settle where they were ultimately observed.
4.3. How good are these models?
The fraction of variation explained by the models (Table 3, right-hand column) ranges
from 10 to 79%, with a mean value of about 33%. In other words, about 67% of the
variation of bivalve counts was not spatially structured or explained by physical
variables or interactions with other size classes. The unexplained variability (21 to 90%)
found in all models suggests phenomena operating at smaller scales, or variability at the
scale of the study site that cannot be expressed as a linear combination of the terms of
the cubic trend-surface equation that we used to model large-scale structures, plus a fair
amount of Poisson sampling error.
Predation by other species, such as waders (Cummings et al., 1997) and eagle rays
(Hines et al., 1997) is a source of variation that we have not examined in this paper. The
name of the embayment where this work has been carried out, Manukau, is a Maori
word meaning ??wading birds?? (X, 1992), by reference to the waders that feed upon
bivalves and other invertebrates. There were consistent patterns of variation in usage by
shorebirds at the scale of sectors (larger than the study site). Effects at the scale of the
study site were expected from previous studies of aggregative responses of shorebirds to
natural variation at the scale of 100 m. Contrary to expectation, there were no increases
in shorebird numbers or large increases in Macomona mortality due to oystercatchers
(Cummings et al., 1997). There was no evidence that shorebirds contributed to the
unexplained variance by substantially altering spatial variation in prey during the
experiment.
It is the extreme size classes (largest and smallest) that were best explained by our
models (last column of Table 3). For the largest animals, the physical variables played
the most important role. For the 0.5?2.5 mm size classes, biological variables (that may
be surrogates for other small-scale physical effects) explained 30 to 40% of the variation
accounted for by the models. The patterns displayed by the smaller animals are both
different from those of the larger animals, and much less spatially structured (Table 2
and Figs. 4?7).
4.4. The physical variables
Elevation and surrogate variables of hydrodynamic forces played a dominant role in
our models. These variables may influence larval deposition, the subsequent transport of
juveniles, and food supply. Physical factors together explained from 10% to about 70%
of the variation found in the bivalve variables across the Wiroa Island study site.
Physical variables derived from hydrodynamic modelling are surrogates for sediment
transport, which is responsible for much of the post-larval dispersal. Grant et al. (1997)
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 123
clearly demonstrate that sediment transport is spatially and temporally variable.
Sediment transport occurs under non-averaged conditions (e.g. high winds) and results
from non-linear interaction between waves, tides, and topography.
In nine bivalve models, elevation was significant, sometimes with its square and cubic
forms (Table 4); many of the bivalve variables showed significant relationships to a
cubic polynomial function of elevation. In most cases, there was a positive partial
2
relationship with elevation, often with a negative relationship with elevation and
3positive with elevation . Elevation controls the action of many of the hydrodynamic
variables that, in turn, may be important determinants of the distributions of the various
bivalve size classes. The next most important variable, significant in eight models, was
the percent of time the plot is covered by more than 20 cm water during spring tide.
Even when the effects of the other variables, including elevation, had been controlled
for, greater numbers of Macomona (all size classes) and Austrovenus (0.5?2.5 mm) were
found in areas covered by water for longer (Table 4). Large Austrovenus (.10 mm)
followed the opposite trend.
The ebb and flood shear stress variables measure the transfer of momentum at the
sediment interface under peak ebb- and flood-tide velocities during a mean tide. Flood
stress made a significant contribution to five models, while ebb stress was significant
only in two. Interestingly, Grant et al. (1997) also found more relationships between
sediment transport and flood stress than ebb stress. Most size classes of both species
were found in greater abundance where flood stress was higher. Small animals tended to
avoid locations with higher values of shell hash (Table 4). The mechanism involved is
unclear since shell hash is buried below the level where juveniles are found. The
distribution of shell hash may be caused by the hydrodynamic processes that also
influence the distribution of juveniles.
The SW and WSW wind-work variables were surrogates for work done on the beach
by wind-driven waves coming from the predominant wind directions. These variables
were only significant for the smallest Austrovenus (22 January), where they may be
surrogates for other unmeasured variables related to elevation. On 22 January, winds
were moderate from the NE, but a 15-knot SW wind had occurred on 19 January. On 16
February, there was a NE breeze up to 15 knots and nothing from the SW. Thus the
modelled wind-wave variables do not perfectly reflect the recent history of the site. On
any given date, recruits may be deposited on the sediment in a pattern determined by
hydrodynamics; changes in the spatial distributions from one time period to another will
depend on the species involved (Hewitt et al., in press).
Adults of Austrovenus, a facultative suspension feeder, occur in the upper 3 cm of
sediment and should be more readily affected by water motion and by physical factors
than large Macomona which live 7?10 cm below the sediment surface. This is precisely
2
what our regression models have shown (larger R for large Austrovenus in step 2 of
Table 3).
4.5. Stability of the bivalve distribution patterns
In agreement with the spatial models, principal component analysis suggests that
124 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
small Macomona and Austrovenus (,4 mm) are controlled by similar processes, since
they are found in the same region of the plane of the first two principal components.
Large Macomona and Austrovenus (.10 mm) are found in opposite locations,
indicating different distributions. Large individuals of each species clearly have the same
distribution on both sampling dates, however (Fig. 8(b, c)). The distributions of
intermediate-sized animals changed from 22 January to 16 February; Macomona (4?15
22
mm) retained the same mean density over the study site (26 animal?m ) but changed
positions (Fig. 8); conversely, the density of Austrovenus (4?10 mm) changed
22dramatically from 2 to 7 animal?m . These findings support the idea that physical
processes are controlling the distribution of this size class.
Larger Macomona and Austrovenus are likely to have lived in the same vicinity for
several years, given their restricted mobility as adults. Regardless of how a cohort of
bivalves ended up inhabiting a given area, if they have been there for several years, their
spatial distribution is a function (likely non-linear) of advection, predation, competition,
etc., i.e. all the factors that could impact them over several years. Historical effects make
it hard to dissect causative factors post hoc, especially compared to smaller bivalves
whose spatial distributions are more likely to bear a close relationship with recent
conditions. If those big bivalves arrived at a smaller size (quite likely, given reduced
passive movement with size) then they will have gone through the same processes as the
smaller ones, but a few years earlier. Thus we are comparing groups of bivalves with
different histories. For that reason, the variables that are likely to determine the spatial
distribution of larger bivalves in our models are those that remain constant or that
represent an integration of physical and biological processes through time, while smaller
bivalves are likely to respond to contemporaneous variables. Elevation, with its second
and third powers, is the variable that dominates the models for larger animals; it can be
seen as an integrated summary of hydrological events over several years.
The potential controlling variables used in this analysis were rarely linearly in-
dependent of one another, as is to be expected in any ecological system, and much of the
initial effort in developing the spatial models involved variable selection. For the
hydrodynamic variables it was necessary to use the results from modelling exercises;
this procedure tends to smooth over extreme values that may be ecologically important
and emphasise average conditions. We have shown that the spatial variation in different
size classes of these bivalves may be related to a variety of ecological processes. The
importance of physical or biological processes varied for different-sized animals, across
a study site chosen to represent basically a flat and physically homogeneous environ-
ment. Our results emphasise the mobility of young bivalves. This mobile post-settlement
phase corresponds to sizes,4 mm for Macomona, and probably even smaller for
Austrovenus (Cummings et al., 1995; Commito et al., 1995b). As they grow, their spatial
distributions slowly modify to adult patterns. However, in contrast to companion studies
(Hewitt et al., 1997; Turner et al., 1997; Thrush et al., 1997a), we found no indication of
adult-juvenile interactions at the scale of our study site. The use of averaged variables,
which may not reflect recent site history, presents difficulties, particularly for the more
mobile juvenile bivalves. Larger individuals are less mobile and integrate large-scale
extrinsic environmental features over long time periods and are thus more suited to this
type of analysis. Variability not explained by the spatial modelling suggests phenomena
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 125
operating at different scales, and the limitation of a linear combination of the terms of
the trend-surface equations, notwithstanding the influence of random factors, local
historical events, and processes operating on different temporal or spatial scales, that
may all blur the picture.
Acknowledgements
We are grateful to the other workshop participants, spouses and friends, who took part
in the Wiroa grid study effort. Special thanks are due to Robert Bell (NIWA) and Roy
Walters (US Geological Survey) who computed data files for several of the physical
variables used in this paper, and checked their descriptions; and to M. J. Anderson
(University of Sydney, NSW, Australia), R. N. Zajac (University of New Haven,
Connecticut, USA) and K. R. Clarke (Plymouth Marine Laboratory, UK) for extensive
comments on preliminary versions of the paper. We thank Auckland Airport Security for
access to the study site. This research was made possible by support from NIWA-NSOF
and FRST-CO1517.
Appendix A
Two examples of spatial modelling (five steps) are presented to illustrate the process.
The first example goes up to step 5 while the second stops at step 4 for lack of
significance. Modelling step 1, using the spatial polynomial equation only, is reported in
Table 2.
1. Macomona 0.5?2.5 mm, 22 January, n5199
2
?Step 2: Physical variables alone: R 50.2234, p,0.0001. y52108.86020.069
Shell hash123.093 Elevation10.044 Ebb stress10.688 Flood stress1107.187
??.20 cm water??
2Step 3: Adding the spatial to the physical variables: no term significant. R for 9
2terms of spatial polynomial50.24662; DR 50.02325.
Step 4: Physical and species variables: attempt to incorporate all the other bivalve
2 2
counts of 22.01: R 50.3785, p,0.0001; DR 50.1551, partial F514.63917,
p,0.0001. 0.406 ??Macomona 2.5?4 mm??10.237 ??Austrovenus 2.5?4 mm??1
0.441 ??Austrovenus 0.5?2.5 mm??
Step 5: Adding the spatial to the physical and biological variables: no term
2 2
significant. R for 9 terms of spatial polynomial50.41045; DR 50.03199.
2. Austrovenus.10 mm, 16 February, n531
2
?Step 2: Physical variables alone: R 50.6872, p,0.0001. y57283.65123189.735
2 3Elevation1856.812 Elevation 2122.907 Elevation 26321.219 ??.20 cm water??
2Step 3: Adding the spatial to the physical variables: R 50.7907, p,0.0001;
126 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
2DR 50.1036, partial F55.93721, p50.0080. Spatial variables (X, Y geographic
2
coordinates) added to the model: 40.725 X 271.248 XY
Step 4: Physical and species variables: attempt to incorporate ??Macomona.15
mm?? of 16.02 in the equation: not significant.
References
Ahn, I.-Y., Lopez, G., Malouf, R., 1993. Effects of the gem clam Gemma gemma on the early post-settlement
emigration, growth and survival of the hard clam Merceneria merceneria. Mar. Ecol. Prog. Ser. 99, 61?70.
Allen, T.F.H., Starr, T.B., 1982. Hierarchy?Perspectives for Ecological Complexity. University of Chicago
Press, Chicago, Illinois, USA, 310 p.
?
? ?Amanieu, M., Legendre, P., Troussellier, M., Frisoni, G.-F., 1989. Le programme Ecothau: theorie ecologique
?et base de la modelisation. Oceanologica Acta 12, 189?199.
Bell, R.G., Hume, T.M., Dolphin, T.J., Green, M.O., Walters, R.A., 1997. Characterisation of physical
environmental factors on an intertidal sandflat, Manukau Harbour. New Zealand. J. Exp. Mar. Biol. Ecol.
216 (1?2), 11?31.
Borcard, D., Legendre, P., Drapeau, P., 1992. Partialling out the spatial component of ecological variation.
Ecology 73, 1045?1055.
Borcard, D., Legendre, P., 1994. Environmental control and spatial structure in ecological communities: an
example using Oribatid mites (Acari Oribatei). Envir. Ecol. Stat. 1, 37?53.
Burrough, P.A., 1987. Spatial aspects of ecological data. In: Jongman, R.H.G., ter Braak, C.J.F., van Tongeren,
O.F.R. (Eds.), Data Analysis in Community and Landscape Ecology. Centre for Agricultural Publishing and
Documentation (Pudoc), Wageningen, The Netherlands, pp. 213?251.
Butman, C.A., 1987. Larval settlement of soft-sediment invertebrates: The spatial scales of pattern explained
by active habitat selection and the emerging role of hydrological processes. Oceanogr. Mar. Biol. Ann. Rev.
25, 113?165.
Cliff, A.D., Ord, J.K., 1981. Spatial Processes: Models and Applications. Pion Limited, London, 266 pp.
Commito, J.A., Currier, C.A., Kane, L.R., Reinsel, K.A., Ulm, I.M., 1995. Dispersal dynamics of the bivalve
Gemma gemma in a patchy environment. Ecol. Monogr. 65, 1?20.
Commito, J.A., Thrush, S.F., Pridmore, R.D., Hewitt, J.E., Cummings, V.J., 1995. Dispersal dynamics in a
wind-driven benthic system. Limnol. Oceanogr. 40, 1513?1518.
Cummings, V.J., Pridmore, R.D., Thrush, S.F., Hewitt, J.E., 1993. Emergence and floating behaviours of
post-settlement juveniles of Macomona liliana (Bivalvia: Tellinacea). Mar. Behav. Physiol. 24, 25?32.
Cummings, V.J., Pridmore, R.D., Thrush, S.F., Hewitt, J.E., 1995. Post-settlement movement by intertidal
benthic macroinvertebrates: do common New Zealand species drift in the water column? N.Z. J. Mar.
Freshwat. Res. 29, 59?67.
Cummings, V.J., Schneider, D.C., Wilkinson, M.R., 1997. Multiscale experimental analysis of aggregative
responses of mobile predators to infaunal prey. J. Exp. Mar. Biol. Ecol. 216 (1?2), 211?227.
Dutilleul, P., Legendre, P., 1993. Spatial heterogeneity against heteroscedasticity: an ecological paradigm
versus a statistical concept. Oikos 66, 152?171.
Emerson, C.W., Grant, J., 1992. The control of soft-shell clam (Mya arenaria) recruitment on intertidal
sandflats by bedload sediment transport. Limnol. Oceanogr. 36, 1288?1300.
Ertman, S.C., Jumars, P.A., 1988. Effects of bivalve siphonal currents on the settlement of inert particles and
larvae. J. Mar. Res. 46, 797?813.
Grant, J., Turner, S.J., Legendre, P., Hume, T.M., Bell, R.G., 1997. Patterns of sediment reworking and
transport over small spatial scales on an intertidal sandflat, Manukau Harbour. New Zealand. J. Exp. Mar.
Biol. Ecol. 216 (1?2), 33?50.
Grant, W.D., Madsen, O.S., 1979. Combined wave and current interaction with a rough bottom. J. Geophys.
Res. 84, 1797?1808.
Haining, R., 1990. Spatial Data Analysis in the Social and Environmental Sciences. Cambridge University
Press, Cambridge, xxi1409 pp.
P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128 127
Haury, L.R., McGowan, J.A., Wiebe, P.H., 1978. Patterns and processes in the time-space scales of plankton
distributions. In: Steele, J.H. (Eds.), Spatial Pattern in Plankton Communities. Plenum Press, New York, pp.
277?327.
He, F., Legendre, P., Bellehumeur, C., LaFrankie, J.V., 1994. Diversity pattern and spatial scale: a study of a
tropical rain forest of Malaysia. Envir. Ecol. Stat. 1, 265?286.
Hewitt, J.E., Pridmore, R.D., Thrush, S.F., Cummings, V.J., in press. Assessing the short-term variability of
spatial patterns of macrobenthos in a dynamic estuarine system. Limnol. Oceanogr.
Hewitt, J.E., Legendre, P., McArdle, B.H., Thrush, S.F., Bellehumeur, C., Lawrie, S.M., 1997. Identifying
relationships between adult and juvenile bivalves at different spatial scales. J. Exp. Mar. Biol. Ecol. 216
(1?2), 77?98.
Hines, A.H., Whitlatch, R.B., Thrush, S.F., Hewitt, J.E., Cummings, V.J., Dayton, P.K., Legendre, P., 1997.
Nonlinear foraging response of a large marine predator to benthic prey: eagle ray pits and bivalves in a New
Zealand sandflat. J. Exp. Mar. Biol. Ecol. 216 (1?2), 191?210.
Isaaks, E.H., Srivastava, R.M., 1989. Applied Geostatistics. Oxford University Press, New York, xix1561 pp.
Legendre, L., Legendre, P., 1983. Numerical Ecology. Elsevier Scientific Publ. Co., Amsterdam, The
Netherlands, xvi1419 pp.
Legendre, P., 1990. Quantitative methods and biogeographic analysis. In: Garbary, D.J., South, R.G. (Eds.),
Evolutionary Biogeography of the Marine Algae of the North Atlantic. NATO ASI Series, G 22,
Springer?Verlag, Berlin, pp. 9?34.
Legendre, P., Fortin, M.-J., 1989. Spatial pattern and ecological analysis. Vegetatio 80, 107?138.
Legendre, P., Borcard, D., 1994. Rejoinder. Envir. Ecol. Stat. 1, 57?61.
Legendre, P., Legendre, L., 1997. Numerical Ecology, 2nd English ed. Elsevier, Oxford.
Legendre, P., McArdle, B.H., 1997. Comparison of surfaces. Oceanologica Acta 20, 27?41.
Luckenbach, M.W., 1984. Settlement and early post-settlement survival in the recruitment of Mulinia lateralis
(Bivalvia). Mar. Ecol. Prog. Ser. 17, 245?250.
Monismith, S.G., Koseff, J.R., Thompson, J.K., O?Riordan, C.A., Nepf, H.M., 1990. A study of model bivalve
siphonal currents. Limnol. Oceanogr. 35, 680?696.
Oden, N.L., 1984. Assessing the significance of a spatial correlogram. Geographical Analysis 16, 1?16.
O?Riordan, C.A., Monismith, S.G., Koseff, J.R., 1993. A study of concentration boundary-layer formation over
a bed of model bivalves. Limnol. Oceanogr. 38, 1712?1729.
Ripley, B.D, 1981. Spatial Statistics. John Wiley and Sons, New York, x1252 pp.
Roegner, C., Andre, A., Lindegarth, M., Eckman, J.E., Grant, J., 1995. Transport of recently settled soft-shell
clams (Mya arenaria L.) in laboratory flume flow. J. Exp. Mar. Biol. Ecol. 187, 13?26.
Schneider, D.C., 1987. Patch structure of benthic populations on an intertidal sandflat. Oceanologica Acta 10,
469?473.
Schneider, D.C., 1994. Quantitative Ecology: Spatial and Temoral Scaling. Academic Press, San Diego,
xv1395 pp.
Snelgrove, P.V.R., 1994. Hydrodynamic enhancement of invertebrate larval settlement in microdepositional
environments: colonization tray experiments in a muddy habitat. J. Exp. Mar. Biol. Ecol. 176, 149?166.
Snelgrove, P.V.R., Butman, C.A., 1994. Animal?sediment relationships revisited: cause versus effect.
Oceanogr. Mar. Biol. Ann. Rev. 32, 111?177.
Sokal, R.R., Oden, N.L., 1978. Spatial autocorrelation in biology. 1. Methodology. Biol. J. Linnean Soc. 10,
199?228.
Stommel, H., 1963. Varieties of oceanographic experience. Science 139, 572?576.
Student [W.S. Gosset], 1914. The elimination of spurious correlation due to position in time and space,
Biometrika, 10, 179?180.
Thrush, S.F., Pridmore, R.D., Hewitt, J.E., Cummings, V.J., 1992. Adult infauna as facilitators of colonization
on intertidal sandflats. J. Exp. Mar. Biol. Ecol. 159, 253?265.
Thrush, S.F., Pridmore, R.D., Hewitt, J.E., Cummings, V.J., 1994. The importance of predators on a sandflat:
interplay between seasonal changes in prey densities and predator effects. Mar. Ecol. Prog. Ser. 107,
211?222.
Thrush, S.F., Hewitt, J.E., Pridmore, R.D., Cummings, V.J., 1996. Adult / juvenile interactions of infaunal
bivalves: contrasting outcomes in different habitats. Mar. Ecol. Prog. Ser. 132, 83?92.
128 P. Legendre et al. / J. Exp. Mar. Biol. Ecol. 216 (1997) 99 ?128
Thrush, S.F., Pridmore, R.D., Bell, R.G., Cummings,V.J., Dayton, P.K., Ford, R., Grant, J., Hewitt, J.E., Hines,
A.H., Hume, T.M., Lawrie, S.M., Legendre, P., McArdle, B.H., Morrisey, D., Schneider, D.C., Turner, S.J.,
Walters, R., Whitlatch, R.B., Wilkinson, M.R., 1997a. The sandflat habitat: scaling from experiments to
conclusions. J. Exp. Mar. Biol. Ecol. 216 (1?2), 1?9.
Thrush, S.F., Cummings, V.J., Dayton, P.K., Ford, R., Grant, J., Hewitt, J.E., Hines, A.H., Lawrie, S.M.,
Legendre, P., McArdle, B.H., Pridmore, R.D., Schneider, D.C., Turner, S.J., Whitlatch, R.B., Wilkinson,
M.R., 1997b. Matching the outcome of small-scale density manipulation experiments with larger scale
patterns: an example of bivalve adult / juvenile interactions. J. Exp. Mar. Biol. Ecol. 216 (1?2), 153?169.
Turner, S.J., Grant, J., Pridmore, R.D., Hewitt, J.E., Wilkinson, M.R., Hume, T.M., Morrisey, D.J., 1997.
Bedload and water-column transport and colonization processes by mobile post-settlement benthic
macrofauna: Does infaunal density matter? J. Exp. Mar. Biol. Ecol. 216 (1?2), 51?76.
Upton, G., Fingleton, B., 1985. Spatial Data Analysis by Example. Vol. 1: Point Pattern and Quantitative data.
John Wiley and Sons, Chichester, xi1409 p.
X, 1992. Maori Place Names, Their Origins and Meanings. Reed Books, Auckland, 95 pp.