Forgot your password?
Document Actions
Anderson and Millar 2004
Journal of Experimental Marine Biology and Ecology
305 (2004) 191 – 221
www.elsevier.com/locate/jembe
Spatial variation and effects of habitat on temperate
reef fish assemblages in northeastern New Zealand
Marti J. Anderson *, Russell B. Millar
Department of Statistics, Tamaki Campus, University of Auckland, Private Bag 92019, Auckland, New Zealand
Received 15 January 2003; received in revised form 5 November 2003; accepted 22 December 2003
Abstract
Reef-associated fishes can respond to changes in habitat structure and the nature of their
response can change with different spatial scales of observation. A structured hierarchical
mensurative sampling design was used to sample temperate reef fish assemblages in northeastern
New Zealand at several spatial scales over 2 years. The three spatial scales examined were tens
of meters (transects), hundreds to thousands of meters (sites) and hundreds of kilometers
(locations). We tested the hypothesis that fish assemblages differed between kelp forest habitat
(relatively dense stands of the kelp, Ecklonia radiata (C. Agardh) J. Agardh, median depth=13.5
m) and barrens habitat (rocky reef dominated by turfing and encrusting red algae and the
grazing urchin, Evechinus chloroticus (Valenciennes), median depth=6.7 m). Recently developed
multivariate techniques were used to test for and quantify multivariate variation at different
spatial scales. There were significant effects of habitat on the spatial distribution of fish
assemblages, characterised by greater abundances or frequencies of Parika scaber, Chromis
dispilus, Trachurus novaezelandiae, Nemadactylus douglasii, Bodianus unimaculatus, Odax
pullus and Pseudolabrus miles in kelp forest habitat, and greater abundances or frequencies of
Notolabrus celidotus, Notolabrus fucicola, Girella tricuspidata, Coris sandageri, Chironemus
marmoratus, Parma alboscapularis, Scorpis violaceus and Kyphosus sydneyanus in barrens
habitat. Some of the more common species, including Upeneichthys lineatus, Scorpis lineolatus
and Cheilodactylus spectabilis showed no strong consistent effects of these two differing habitats
on their distributions. There was, however, a significant HabitatÂLocations interaction: effects of
habitat did not occur at all locations. Variability was highest at the scale of individual transects
and variability from site to site and from location to location was comparable. Spatial variation
was large compared to inter-annual variation, which was minimal, and spatial patterns were
consistent in the 2 years examined. Further experiments, including manipulations, are required to
* Corresponding author. Tel.: +64-9-373-7599x85052; fax: +64-9-373-7000.
E-mail address: mja@stat.auckland.ac.nz (M.J. Anderson).
0022-0981/$ - see front matter D 2004 Elsevier B.V. All rights reserved.
doi:10.1016/j.jembe.2003.12.011
192 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
understand what mechanisms and processes might be driving these patterns. This study, coupled
with results from previous studies, suggests that there may be a dynamic inter-play between
effects of habitat on fish and effects of fish on biogenic habitat, such as kelp forests.
D 2004 Elsevier B.V. All rights reserved.
Keywords: Hierarchical experimental design; Kelp forest; Multivariate analysis; Temperate reef fish; Urchin
grazing; Variance components
1. Introduction
An important goal in ecology is to understand patterns of distributions of organisms
by reference to the habitat available to them in the environment. For reef-associated
fishes, there is abundant evidence that the structure of the habitat has important effects
on spatial distributions of populations, both in tropical coral reefs (Roberts and
Ormond, 1987; Tolimieri, 1995; Caley and St. John, 1996; Friedlander and Parrish,
1998; Tolimieri, 1998a; Holbrook et al., 2000; McClanahan and Arthur, 2001) and in
temperate rocky reef systems (Choat and Ayling, 1987; Jones, 1988; Holbrook et al.,
1990; Connell and Jones, 1991; Carr, 1989; Levin and Hay, 1996; Tupper and
´ ´
Boutilier, 1997; Garcıa-Charton and Perez-Ruzafa, 2001). For example, for temperate
reefs in northeastern New Zealand, Choat and Ayling (1987) described two distinct
assemblages of fishes associated with two different habitat types: (i) areas dominated
by grazing urchins (Evechinus chloroticus (Valenciennes)), with cover by encrusting
and turfing red algae (called ‘‘barrens’’ habitat) and (ii) areas of relatively dense
stands of the laminarian kelp Ecklonia radiata (C. Agardh) J. Agardh (referred to as
‘‘kelp forests’’). They found that small wrasses were more abundant in the shelter of
kelp forests and increased with increases in kelp density, while larger carnivores were
more abundant in urchin-grazed areas (Choat and Ayling, 1987). Jones (1984a) found
that juveniles of the wrasse, Notolabrus celidotus (spotty), were positively associated
with the density of macroalgae on New Zealand reefs. Furthermore, densities of
recruits decreased with the experimental removal of kelp and increased with the
experimental addition of kelp (through removal of urchins, Jones, 1984a). The
consequences of the use of kelp forests by younger fish for predicting distributions
of adult populations are, however, largely unknown.
There is some evidence for the effects of kelp habitat on reef fishes in other parts
of the world. However, Holbrook et al. (1990) found only weak differences in species
composition among reefs with different types of algal habitats in Southern California,
and mixed results were obtained by other workers in this region (e.g., Stephens et al.,
1984; Carr, 1989; DeMartini and Roberts, 1990). In addition, in New Zealand, Choat
and Ayling (1987) suggested that the differences found between habitats in their study
were independent of species identity, being driven instead by differences in the
biology of feeding preferences of fishes at different ontogenetic stages. Thus, it is
unclear whether the species composition and relative abundances of individual fish
species differ in consistent or predictable ways between these two identifiable habitats
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 193
in northeastern New Zealand. A structured mensurative experiment is needed to
investigate this.
A further important aspect of understanding spatial distributions with respect to
habitat characteristics is to recognize that observed patterns are dependent on the
spatial scale of observation (e.g., Andrew and Mapstone, 1987; Wiens, 1989; Tolimieri,
1995; Chesson, 1998; Sale, 1998). There have been several studies of spatial patterns
of distributions of fish at several spatial scales, from meters up to hundreds or
thousands of kilometers (e.g., Choat and Ayling, 1987; Doherty, 1987; Fowler et al.,
1992; Tolimieri, 1998a; Ault and Johnson, 1998; Connell and Kingsford, 1998; Garcıa- ´
´
Charton and Perez-Ruzafa, 2001; Connell, 2002). Effects of habitat at small spatial
scales can provide predictive power at larger spatial scales for some coral reef fish
species (i.e., damselfish, Tolimieri, 1995; Holbrook et al., 2000), but not for others
(i.e., stoplight parrotfish, Tolimieri, 1998a,b). Jones (1988) suggested that it is the
comparison of the magnitude of variation at each spatial scale of interest that will
provide the necessary framework for formulating hypotheses about relevant processes
for reef fishes. Hierarchical spatially structured sampling programs provide a means of
partitioning and quantifying the magnitude of variation at different spatial scales
(Andrew and Mapstone, 1987; Underwood and Chapman, 1996; Underwood et al.,
2000). Indeed, it may not be possible to understand effects of habitat, which is a
spatial phenomenon, involving patchiness and heterogeneity, without an understanding
of variability at different spatial scales (Kotliar and Wiens, 1990; Underwood and
Chapman, 1996; Underwood et al., 2000). For example, Fowler-Walker and Connell
(2002) demonstrated scale-dependent associations of understorey algae in E. radiata
forests, with weak or variable local-scale responses to habitat, but strong and consistent
patterns at regional scales (thousands of kilometers). Although some explicit measure-
ments of variation in effects of habitat on fish species at different spatial scales have
been made by certain workers (e.g., Tolimieri, 1998a; Doherty, 1987; Fowler et al.,
´ ´
1992; Garcıa-Charton and Perez-Ruzafa, 2001), this has not generally been done using
structured hierarchical designs, such as those used to advantage for benthic invertebrate
and algal assemblages (e.g., Underwood and Chapman, 1996; Menconi et al., 1999;
Benedetti-Cecchi, 2001; Fowler-Walker and Connell, 2002), but see Connell and
Kingsford (1998) and Connell (2002).
One stumbling block towards measuring spatial variation in fish or other kinds of
assemblages has been the lack of an appropriate method for assessing multivariate
variation for several species simultaneously, that is, to measure and quantify variation
in multivariate assemblages in each scale of a hierarchy of spatial scales. Available
distance-based multivariate methods that are realistic for non-normal counts of species
abundances are either unable to obtain independent partitions of multivariate variation for
such complex designs (such as ANOSIM, Clarke, 1993; or Mantel correlograms, Legendre
and Fortin, 1989), or require large numbers of replicates to avoid problems of non-
independence (Underwood and Chapman, 1998).
Recent developments in non-parametric multivariate analysis provide a method for
analyzing multivariate assemblages on the basis of any distance or dissimilarity
measure, while also allowing the data to be partitioned according to any experimental
design, including nested hierarchies (Anderson, 2001a; McArdle and Anderson, 2001).
194 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Non-parametric multivariate analysis of variance (NPMANOVA, Anderson, 2001a) can
be used to partition variability and to provide measures of multivariate variability at
different scales in a structured hierarchical design, in a manner directly analogous to
univariate partitioning using ANOVA (e.g., Searle et al., 1992; Benedetti-Cecchi,
2001). The statistical significance of multivariate variance components can also be
tested using appropriate methods of permutation (Anderson, 2001b; Anderson and ter
Braak, 2003).
The purpose of the present investigation was to conduct a mensurative, observa-
tional experiment (sensu Hurlbert, 1984) to examine the potential effects of habitat on
fish assemblages in northeastern New Zealand at several spatial scales. Before
hypotheses concerning any potential emerging patterns or processes can be developed,
good observational data and quantitative measures of spatial variability in fish
assemblages are extremely useful—for individual species and for the multivariate
assemblage as a whole (Underwood et al., 2000). More specifically, we wished to test
the null hypothesis (following Choat and Ayling, 1987) that there are no differences in
the composition of fish assemblages in kelp forest habitat vs. barrens habitat on
subtidal rocky reefs (i.e., that differences are ‘‘independent of species identity’’). To
determine the generality of any potential patterns, effects of habitat need to be
examined at several locations and at several times (e.g., Yates and Cochran, 1938;
Snedecor, 1946; Underwood and Petraitis, 1993). Furthermore, as effects of habitat on
fish may or may not ‘‘scale up’’ to provide good prediction at larger spatial scales
(e.g., Tolimieri, 1998a), depending on the species (e.g., Ault and Johnson, 1998), we
predicted that the multivariate effects of habitat would interact with variability at
different spatial scales and that different species would show different patterns in this
regard. We also tested the hypothesis that there is a significant relationship between
fish assemblage structure and the density of kelp forests, as a relationship with kelp
density was observed for total numbers of large carnivorous fish by Choat and Ayling
(1987).
The spatial scales of investigation included in the study were tens of meters
(transects), hundreds of meters to kilometers (sites) and hundreds of kilometers
(locations). These large-scale surveys were performed in each of 2 years using a
structured hierarchical experimental design to investigate effects of habitat on reef
fishes. This is the first study, to our knowledge, that uses multivariate methods to obtain
independent quantitative measures of variability in fish assemblages at different spatial
scales.
In northeastern New Zealand, it is known that patches of kelp forest habitat tend to
occur at deeper depths than patches of barrens habitat, although there is some overlap in
their natural depth ranges (e.g., Schiel, 1990). Brook (2002) found a significant positive
relationship between depth (up to 45 m) and species richness of fishes, while depth-
stratified sampling by Meekan and Choat (1997) showed significant differences in
abundances of several prominent herbivorous fishes among different depth strata. It is
not known, however, the extent to which effects of depth may mediate effects of habitat.
The present study therefore also measured depth as a covariable for analyses, to test the
hypothesis that significant differences in assemblages due to habitat, if present, could not
be fully explained by a relationship of assemblages with depth.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 195
2. Methods
2.1. Underwater visual sampling of fish
The study used a structured hierarchical experimental design. Four locations, separated
by hundreds of kilometers, were sampled along the northeastern coast of New Zealand
(Fig. 1). These included, from South to North, Hahei (36j50.86VS, 175j49.32VE), Leigh
(36j17.43VS, 174j48.82VE), Home Point (35j19.38VS, 174j21.38VE) and Berghan Point
(34j55.78VS, 173j32.72VE) (Fig. 1). Within each location, fish were sampled at each of
four different randomly located sites, separated by hundreds to thousands of meters. At
each site, two habitats were investigated: kelp forests (i.e., areas characterised by relatively
dense cover of the kelp E. radiata) and ‘‘barrens’’ (i.e., areas characterised by little or no
macroalgal cover and dominated by the grazing urchin E. chloroticus). The range of
depths was 2 –20 m, although the majority of the observations were in the range from 5 to
15 m.
Within each habitat, divers on SCUBA did a visual survey by swimming along a
transect and identifying and recording the number of each species of fish observed
within a distance of 2.5 m on either side of the transect. Taxonomic authorities for all
fish species named herein are available in Paulin et al. (1989). All fish seen were
recorded except for notoriously cryptic or very small species (e.g., Tripterygiidae and
Gobiidae), which were not included, as we could not be confident that they were being
reliably seen over these spatial scales (e.g., Lincoln-Smith, 1989; Willis, 2001). Each
diver carried a tape measure and swam a ‘‘run-in’’ distance of 5 m before beginning the
actual survey for each transect, which was 25 m long after the run-in. Thus the total area
sampled per transect was 5Â25 m. There were n=10 transects sampled haphazardly
within each habitat at each site. Although there are some known biases inherent in using
the method of underwater visual surveys to count fish (e.g., Brock, 1982; Sale and
Sharp, 1983; Lincoln-Smith, 1989; Thompson and Mapstone, 1997; Willis et al., 2000),
this same sampling methodology was used across the entire experimental design. This
allows valid comparisons to be made across habitats and across different spatial scales,
even though unbiased estimates of population densities may not be obtained for some
fish species from these data. It was not possible logistically for a single diver to perform
the entire sampling design, so several divers participated in the study. There was likely
variation due to different observers, although this was not explicitly measured. This was
not, however, related in any systematic fashion with differences across habitats because
each diver involved in the study sampled both types of habitat. Divers also recorded the
depth at the beginning and at the end of each transect (in meters), which generally
differed by no more than 1 – 2 m. The average depth measured for each transect was
used for analyses. Also, in the first year of sampling for kelp habitats, the density of
kelp in the forests was estimated. This was done by a second diver swimming behind
the diver counting fish. A 1-m2 quadrat was placed alongside the transect tape (used by
the diver counting fish) at each of five positions along the tape (approximately 5, 10, 15,
20 and 25 m) and the number of kelp plants of the species E. radiata per quadrat was
recorded. The average and the standard deviation for kelp density were then calculated
for each transect based on these five observations.
196 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Fig. 1. Map of northeastern New Zealand showing the four locations for the study and the sites where surveys
were done. Grey circles indicate kelp habitat and black triangles indicate barrens habitat sampled at each site.
The sampling was done within a 1-month period in each of 2 consecutive years during
the southern hemisphere’s (austral) summer: from 30 November to 21 December 2000
(year 1) and from 7 January to 5 February 2002 (year 2). The locations were not sampled
consecutively from north to south or vice versa, so as not to confound the latitudinal
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 197
gradient (if any) with time of sampling within the month. In addition, all sites were located
outside of any established marine reserves (i.e., at Hahei and Leigh). GPS coordinates
were recorded at each site in year 1 and the same sites (although not the same transects)
were re-visited in year 2 by reference to these GPS coordinates. Due to the high mobility
of fish, the 2 years were considered as independent samples.
2.2. Multivariate statistical methods
The experimental design therefore consisted of four factors: year (two levels, fixed),
location (four levels, random), site (four levels, random, nested in location) and habitat
(two levels, fixed, crossed with all other factors). Thus, with n=10 transects, there were a
total of 640 observation units in the data set. There were 46 species of fish (variables)
recorded in the study that were included in multivariate analyses (listed in Appendix A).
Although it would have been logical to treat the factor ‘‘year’’ as random in the present
design, with only 2 years of data, inferences concerning inter-annual variation in general
could not be viewed as very precise. Thus, ‘‘year’’ was treated as a fixed factor, which
focused our analyses on results obtained for these 2 years only. With a greater number of
years of sampling anticipated in the future, we would treat this as a random effect.
The fish species variables were highly skewed and contained a great many zeros,
making traditional analyses (which assume normality of errors) unsuitable, so non-
parametric approaches were used. NPMANOVA (Anderson, 2001a; McArdle and Ander-
son, 2001) was used to analyse the multivariate data set in response to the complete
experimental design (including interactions). This method allows multivariate data to be
analysed on the basis of any distance or dissimilarity measure of choice, in response to any
complex experimental design, with P-values obtained using permutations. We used scale-
invariant binomial deviance as a dissimilarity measure (described below) for the NPMA-
NOVA. The tests do not assume that the original variables conform to a multivariate
normal or to any other specified distribution. Although the approach does not explicitly
assume common variances among groups, it will (like ANOSIM, Clarke, 1993) be
sensitive to differences in multivariate dispersion. In addition, the results of NPMANOVA
were also used to estimate the sizes of multivariate pseudo variance components from the
analyses, relying on the analogy of NPMANOVA as an ANOVA based on distances and
univariate ANOVA estimators based on mean squares (e.g., Searle et al., 1992). Doing this
(and indeed applying NPMANOVA to a complex multi-factorial design in general) does
require us to assume that the linear ANOVA model can be applied successfully to achieve
a partitioning of the squared inter-point dissimilarities for inferences to be made about the
multivariate assemblage. We feel this is reasonable provided the choice of dissimilarity
measure accurately reflects relevant qualities of the multivariate assemblage of interest.
For each term in the analysis, 4999 permutations of raw data units were done to obtain
P-values (e.g., Manly, 1997). Care was taken to ensure that the correct permutable units
were used to obtain a valid permutation test of each term in the analysis (Anderson,
2001b; Anderson and ter Braak, 2003). For example, to test the following term:
YearÂLocationÂHabitat, the 64 cell units corresponding to the 64 combinations of
YearÂSite(Location)ÂHabitat were permuted (i.e., the 10 transects within each of these
cells were kept together as a unit). This is because the denominator mean square for the
198 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
test of YearÂLocationÂHabitat is the mean square for YearÂSite(Location)ÂHabitat,
thus it is these units (indicated by the denominator of the F-ratio) that are appropriate to
permute under the null hypothesis (Anderson, 2001b; Anderson and ter Braak, 2003). An
important assumption underlying each test is therefore the exchangeability of the correct
errors (for each particular term being tested in the analysis) under each null hypothesis of
interest. In some cases, this restriction on the permutable units meant that there were not
enough possible permutations to get a reasonable test. For these situations, a Monte Carlo
sample was drawn from the theoretical asymptotic permutation distribution (Anderson
and Robinson, 2003). We re-iterate that these permutation tests assume not just
exchangeability, but also that a linear model on dissimilarities is appropriate for choosing
reasonable test-statistics for testing multivariate hypotheses in an ANOVA framework,
analogous to those used in univariate ANOVA. All analyses were done using specialised
software, written in FORTRAN.
To compare whole fish assemblages to quantitative variables (i.e., depth, average kelp
density and the standard deviation of kelp density), non-parametric multivariate multiple
regression was used on the basis of the binomial deviance dissimilarity measure, using
4999 random permutations (McArdle and Anderson, 2001). We were particularly
concerned to test for any significant effects of habitat over and above the potential effect
of depth, as kelp forests and barrens habitats occur naturally at different depths, on
average. Thus, the effect of habitat was tested, with depth as a covariable in the analysis,
using non-parametric multivariate multiple regression and 4999 permutations of residuals
under the reduced model (e.g., Anderson, 2001b). We also used ordinary least squares
regression to compare each of the three quantitative variables to the total number of
species per transect and the total number of individuals (the latter after transforming to
ln(x+1)). Individual transects were removed from the analysis whenever the quantitative
predictor variable of interest was missing. Also, for these and any other analyses of the
total number of fish (as a univariate variable), the following schooling species that can
occur sporadically in clumps of thousands were not included: Trachurus novaezelandiae
(jack mackerel), Aldrichetta forsteri (yellow-eyed mullet) and Decapterus koheru
(koheru). Far from being dominant species, mullet and koheru, for example, were only
observed in 4 and 16 out of 640 transects, respectively. More importantly, these species
were excluded from sums of the total number of individuals primarily because no certainty
regarding the precision of estimates of their numbers was warranted (Choat and Ayling,
1987). These species were not, however, omitted from the multivariate analyses, as such
large numbers do not pose a problem for the dissimilarity measure used here. Although
there still may be some difficulty in interpreting analyses of total numbers of individuals
among sites with different species composition, as the underwater visual sampling method
may be biased in different ways for different fish species, we nevertheless considered that
it was a measured variable worth investigating, with this caveat in mind.
2.3. Binomial deviance as a dissimilarity measure
For the fish abundance data, we developed and used a new dissimilarity measure, based
on likelihood theory. Let y1k and y2k be the count for species k in transects 1 and 2,
respectively, and let nk=( y1k+y2k). We may then consider the null hypothesis that the two
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 199
transects do not differ in their composition and relative abundance of species. That is, for
each species, we expect that half of the counts (i.e., half of nk) will fall into each transect.
The binomial deviance of our observed data from this hypothesis is defined, under
likelihood theory, as
y1k y2k 1
deviance ¼ y1k log þ y2k log À ðy1k þ y2k Þlog
nk nk 2
Thus, a useful measure of ecological dissimilarity between the two transects is the sum of
these deviances across all species. To take account of the fact that different species may be
varying on different scales, a scale-invariant measure may be obtained by considering this
quantity on a per-observation basis, which is achieved by simply dividing by nk:
X1
p
y1k
Y2k
1
d1;2 ¼ y1k log þ y2k log À ðy1k þ y2k Þlog
k¼1
nk nk nk 2
This measure, which we will refer to as the binomial deviance dissimilarity, resembles
somewhat the measure proposed by Cao et al. (1997) as an improvement on the Bray-
Curtis measure, but has the added advantage of being based in likelihood theory.
The null hypothesis from which the binomial deviance dissimilarity is derived could
also be tested using the chi-square test statistic for the equality of expected counts in the
two transects for all species. In this context, the chi-square test statistic is equivalent to the
dissimilarity measure known as the coefficient of divergence (Clark, 1952). Hence, the
binomial deviance and the coefficient of divergence should have similar properties. The
coefficient of divergence, along with the Canberra metric and the Bray-Curtis measure,
were recommended by Gower and Legendre (1986) for use with species abundance data.
For further details on the properties of these measures, see Legendre and Legendre, (1998,
e.g., p. 298).
2.4. Ordination
Metric multi-dimensional scaling (MDS) (principal coordinate analysis) on the basis of
the binomial deviance dissimilarity measure was used as an unconstrained ordination
method to visualize multivariate patterns. In addition, canonical analysis of principal
coordinates (CAP, Anderson and Willis, 2003; Anderson and Robinson, 2003) was used as
a constrained ordination procedure to visualize patterns by reference to particular
hypotheses. The CAP analyses were only done on appropriate terms found to be
significant using NPMANOVA. These analyses were also done using specialised software
written in FORTRAN.
2.5. Univariate analyses
Frequencies of occurrence were examined for each fish species across the relevant
factors of interest in the study (Appendix A). Six species of fish were abundant enough to
be analysed using univariate analysis of variance (ANOVA): Parika scaber (leather
200 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
jacket), Cheilodactylus spectabilis (red moki), Scorpis lineolatus (sweep), Upeneichthys
lineatus (goatfish), N. celidotus (spotty) and Chromis dispilus (demoiselle). The total
number of fish and the total number of species per transect were also analysed using
univariate ANOVA. Due to the predominance of zeros, extremely large variability at the
level of individual transects and the highly skewed nature of these data, normality was not
a reasonable assumption for any of these single variables except for the total number of
species. Thus, in each case, all tests were done using a permutation procedure (with 4999
permutations of appropriate units), as described for the NPMANOVA above. These tests
using permutations will be sensitive to differences in dispersion, so results should be
interpreted with caution in this regard. Where appropriate, significant terms were
investigated with a posteriori pair-wise comparisons, which also used 4999 random
permutations to obtain P-values. These univariate tests were done using the same software
as that used for the multivariate tests, but with only one variable and (appropriately) basing
the analysis on Euclidean distance. The F-ratios used for tests done in this way are
equivalent to those of traditional ANOVA (Anderson, 2001a), although the P-values are
not obtained using traditional tables.
3. Results
3.1. Effects of depth and density of kelp
The two habitats sampled did differ in their spatial distributions with respect to depth
(Fig. 2). The median depth of transects in kelp habitat was 13.5 m (inter-quartile
range = 11.5 –15 m), while for barrens habitat the median depth was 6.7 m (inter-quartile
Fig. 2. Boxplots of the depth of transects (in metres) for each habitat: barrens (n=316) and kelp (n=318).
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 201
Table 1
Relationships between the total number of fish species (tot sp) and the total number of fish (tot fish) vs. depth,
average kelp density (kelp av) or standard deviation of kelp density (kelp sd)
n r2 F P
Year 1
tot sp vs. depth 309 0.0043 1.325 0.2506
tot fish vs. depth 309 <0.0001 0.003 0.9538
tot sp vs. kelp av 159 0.0711 12.022 0.0007
tot fish vs. kelp av 159 0.0167 2.700 0.1043
tot sp vs. kelp sd 159 0.0774 13.171 0.0004
tot fish vs. kelp sd 159 0.0349 5.673 0.0184
Year 2
tot sp vs. depth 325 0.0057 1.869 0.1726
tot fish vs. depth 325 0.0157 2.691 0.0238
The total number of fish was transformed to ln(x+1) before analysis. Note that kelp densities were not recorded in
year 2.
range = 5.35– 8 m). There was no significant relationship, however, between the total
number of species and depth for either year 1 or for year 2 (Table 1). The total number of
fish (transformed to ln(x+1)) was significantly and positively related to depth for year 2
only (Table 1). There was also a significant relationship between the depth of transects and
the multivariate fish assemblages (Table 2). Depth only explained 3.5% of the variation in
the multivariate assemblage structure and only 1.6% of the variation in total numbers of
fish. Furthermore, effects of habitat (kelp vs. barrens) were significant over and above
effects of depth (Table 2). Thus, although some effects of habitat (described in greater
detail below) may be attributable to differences in depth, this analysis shows there were
significant effects of habitat on fish assemblages (e.g., structural or other differences) that
were unrelated to depth. Depth was not treated as a covariable in subsequent analyses (e.g.,
Table 3 below), as we wished to examine overall effects of habitat, including those aspects
that may have been related to depth.
Within kelp forest habitats for year 1, fish assemblages were also significantly related to
the average kelp density per transect (pseudo F1, 157=4.788, P=0.0008, 4999 permuta-
tions), but were not related to the variation (standard deviation) in kelp density along
transects (pseudo F1, 157 = 1.214, P = 0.3134, 4999 permutations). The total number of
Table 2
Sequential non-parametric multivariate multiple regression showing the relationship between multivariate fish
species abundance data (based on the binomial deviance dissimilarity measure) and depth, followed by the effect
of habitat, taking depth into account as a covariable
Source df SS MS F P
Depth 1 32.660 32.660 22.942 0.0002
Habitat/depth 1 10.050 10.050 7.128 0.0078
Residual 631 889.665 1.410
Total 633 932.375
202 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 3
Non-parametric multivariate analysis of variance of 46 fish species abundance variables, based on the binomial
deviance dissimilarity measure
Source df SS MS F P Denominator Permutable
MS units
Year (Ye) 1 8.280 8.280 0.979 0.4608 YeÂLo 8 YeÂLo units
Location = (Lo) 3 87.368 29.123 5.765 0.0002 Si(Lo) 16 Si(Lo) units
Site(Location) = 12 60.620 5.052 4.388 0.0002 Res 640 observation
Si(Lo) units
Habitat = Ha 1 25.708 25.708 2.264 0.0650 LoÂHa 8 LoÂHa units
YeÂLo 3 25.384 8.461 2.975 0.0012 YeÂSi(Lo) 32 YeÂSi(Lo)
units
YeÂSi(Lo) 12 34.130 2.844 2.470 0.0002 Res 640 observation
units
YeÂHa 1 4.117 4.117 1.598 0.1948 YeÂLoÂHa 16 YeÂLoÂHa
units
LoÂHa 3 34.065 11.355 3.684 0.0002 Si(Lo)ÂHa 32 Si(Lo)ÂHa
units
Si(Lo)ÂHa 12 36.983 3.082 2.677 0.0002 Res 640 observation
units
YeÂLoÂHa 3 7.728 2.576 1.256 0.2248 YeÂSi(Lo)ÂHa 64 YeÂSi(Lo)Â
Ha units
YeÂSi(Lo)ÂHa 12 24.609 2.051 1.781 0.0004 Res 640 observation
units
Residual 576 663.153 1.151
Total 639 1012.145
P-values were obtained using 4999 permutations of given permutable units for each term or using 4999 Monte
Carlo samples from the asymptotic permutation distribution (given in italics) when there were few possible
permutations.
species was significantly positively related to the average and to the standard deviation in
kelp density per transect, while the total number of fish (transformed to ln(x+1)) was
significantly positively related only to the standard deviation in kelp density (Table 1).
However, these relationships, although statistically significant, were very weak, with a
large amount of scatter and very low values of r2 (Table 1). This indicates that kelp density
and depth, while apparently contributing to explain small amounts of observed variation in
fish assemblages, do not provide the means to generate predictions at the scale of
individual transects.
3.2. Measured variation
Non-parametric MANOVA on the fish data showed that there was significant small-
scale variability in the fish assemblages from site to site and year to year in different
habitats (i.e., a significant YearÂSite(Location)ÂHabitat interaction, Table 3). This
significant small-scale variation was reflected in the analysis of several of the
relatively abundant individual species as well, namely: N. celidotus, S. lineolatus
and C. dispilus (see Tables 6 and 8 below). The greatest multivariate variability
occurred at the scale of individual transects (i.e., residual, Table 4). The next-greatest
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 203
Table 4
Estimated multivariate pseudo variance components for each term in the model based on sums of squared
dissimilarities (binomial deviance) for 46 fish species and the analogous univariate ANOVA estimators using
mean squares (i.e., using multivariate mean squares in Table 1), ordered from largest to smallest
Source MS Variance components
Residual 1.1513 1.1513
Lo 29.1227 0.1504
LoÂHa 11.3551 0.1034
Si(Lo) 5.0516 0.0975
Si(Lo)ÂHa 3.0819 0.0965
YeÂSi(Lo)ÂHa 2.0507 0.0899
YeÂSi(Lo) 2.8441 0.0846
Ha 25.7084 0.0449
YexLo 8.4614 0.0702
YeÂLoÂHa 2.5760 0.0131
YeÂHa 4.1171 0.0096
Ye 8.2796 0.0000
All components are random variance components except for habitat, year and YeÂHa which are sums of squared
fixed effects.
measured source of variation was that across locations, followed by the random
interaction of location with habitat, and then the variation from site to site within
locations. Considerable variation was also contributed by the random interaction of
sites with habitats (Table 4). Other measured sources of variation included interactions
of spatial effects with years; however, overall changes in assemblages from year to
year were not strong compared to the other effects. In fact, the estimated variation due
to years was slightly negative (À0.0006) using the ANOVA estimator and thus was
estimated as zero (Table 4).
3.3. Effects of habitat
The effects of habitat (kelp forest vs. barrens) on fish assemblages varied significantly
across locations (Table 3). In particular, the effects of habitat appeared to be strongest at
Hahei and Home Point. This is shown by a greater separation of the points representing
assemblages in kelp vs. barrens assemblages for these two locations, compared to Berghan
Point and Leigh in the unconstrained ordinations: the two-factor metric MDS plot (Fig. 3)
and in each of the one-factor metric MDS plots (Fig. 4, left-hand side). In addition, the
CAP analysis showed larger canonical correlation coefficients for Hahei and Home Point
than for Leigh or Berghan Point (Table 5, Fig. 4, right-hand side). Furthermore, the
allocation success showed that fish assemblages from the two habitats at either Leigh or
Berghan Point were not as distinct as that seen for the other locations. For example, Leigh
had only 69% allocation success compared to Hahei, which had 94% success (Table 5).
Pair-wise comparisons showed that there were significant differences in assemblages of
fishes between barrens and kelp habitats at either Hahei or Home Point, but not at Berghan
Point or Leigh (Table 5). It is interesting to note that the CAP analysis for Leigh shows a
fairly distinct separation due to habitat, although its allocation success was poor and there
was no significant effect detected by NPMANOVA. This is essentially due to the selection
204 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Fig. 3. Two-factor metric MDS plot of the pooled fish assemblages for each site on the basis of the binomial
deviance dissimilarity measure, showing the factors of location and habitat. There are eight points for each
combination of the factors, which correspond to the four sites in each of 2 years.
effects occurring in high-dimensional space, which can cause canonical plots to paint a
‘‘rosy’’ picture of the results. It further emphasizes the importance of using a technique
like the leave-one-out method (e.g., Lachenbruch and Mickey, 1968) for assessing the
distinctness of groups and independent tests rather than relying solely on the canonical
plot.
There were no consistent effects of habitat on the total number of fish; these effects
varied from site to site and from year to year (Table 6a, Fig. 5a). Despite the
significant three-way interaction of YearÂLocationÂHabitat (Table 6b), there were no
statistically significant pair-wise differences detected in the total number of species in
kelp forest vs. barrens in either year at any location (pair-wise comparisons, P > 0.05,
Fig. 5b). In contrast, average numbers of P. scaber (leather jackets) were significantly
higher in kelp forests than in barrens habitats for both years at the two northern
locations of Berghan Point and Home Point (Table 6c, Fig. 5c, P < 0.05). Their
frequency of occurrence was also greater in kelp forests, as was that of C. dispilus
(demoiselle), T. novaezelandiae (jack mackerel), Nemadactylus douglasii (porae),
Bodianus unimaculatus (pigfish), Odax pullus (butterfish) and Pseudolabrus miles
Fig. 4. Unconstrained metric MDS plots (left) and constrained CAP plots (right) done separately at each location
(rows) on the basis of the binomial deviance dissimilarity measure, in each case comparing fish assemblages in
two different habitats: barrens vs. kelp. There are eight points for each combination of the factors, which
correspond to the four sites in each of 2 years.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 205
206 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 5
Results of CAP analyses examining effects of habitat within each location for 46 species of fish on the basis of the
binomial deviance dissimilarity measure
Location m %Var Allocation success (%) d2 P
Barrens Kelp Total
Berghan Point 4 78.56 75 75 75 0.607 0.0564
Home Point 5 84.81 88 75 81 0.748 0.0018
Leigh 6 86.97 63 75 69 0.690 0.1080
Hahei 4 86.89 88 100 94 0.824 0.0006
m = the number of principal coordinate (PCO) axes used in the CAP procedure, %Var = the percentage of the total
variance explained by the first m PCO axes, Allocation success = the percentage of points correctly allocated into
each group, d2 = the squared canonical correlation. P-values given are the results of pair-wise comparisons of
assemblages in kelp forest vs. barrens habitat at each location, using NPMANOVA, with 4999 permutations of
individual sites as units (i.e., permuting each set of n = 10 transects together as a unit).
(scarlet wrasse) (Appendix A). Correlations of species with canonical axes also
suggested that Zeus faber (john dory) and Seriola lalandi (kingfish) were associated
with kelp forest habitats (Table 7). Although average numbers of C. dispilus also
tended to be greater in kelp forest habitat (Fig. 5d), this pattern was not statistically
significant (Table 6d).
Fish that occurred more frequently in barrens habitats were N. celidotus (spotty),
Notolabrus fucicola (banded wrasse), Girella tricuspidata (parore), Coris sandageri
Table 6
Results of univariate ANOVAs on selected variables and species
Source df (a) Total no. fish (b) Total no. species (c) Parika scaber (d) Chromis dispilus
MS F P MS F P MS F P MS F P
Ye 1 8917.0 4.779 0.1114 3.025 0.070 0.8108 2.336 0.892 0.4206 4944.4 4.565 0.1164
Lo 3 13056.2 4.065 0.0370 98.173 5.342 0.0170 15.135 2.789 0.0888 9972.5 7.064 0.0088
Si(Lo) 12 3211.7 4.916 0.0002 18.377 4.853 0.0002 5.427 2.674 0.0014 1411.7 4.575 0.0002
Ha 1 17.2 0.010 0.9332 32.400 0.388 0.5752 43.403 2.908 0.1788 3048.9 4.761 0.1084
YeÂLo 3 1866.0 1.934 0.1824 43.050 9.156 0.0026 2.618 1.329 0.3076 1083.2 1.516 0.2600
YeÂSi(Lo) 12 965.1 1.477 0.1268 4.702 1.242 0.2538 1.970 0.970 0.4858 714.4 2.315 0.0068
YeÂHa 1 1196.5 0.586 0.5122 12.656 0.762 0.4538 9.344 8.959 0.0572 3032.5 5.714 0.0948
LoÂHa 3 1794.5 0.564 0.6510 83.508 7.402 0.0032 14.924 6.991 0.0036 640.4 0.738 0.5406
Si(Lo)ÂHa 12 3184.4 4.874 0.0002 11.281 2.979 0.0010 2.135 1.052 0.4020 867.2 2.810 0.0022
YeÂLoÂHa 3 2043.5 1.335 0.3046 16.606 4.516 0.0278 1.043 0.900 0.4670 530.7 0.523 0.6790
YeÂSi(Lo)Â 12 1530.7 2.343 0.0062 3.677 0.971 0.4762 1.159 0.571 0.8830 1013.9 3.286 0.0002
Ha
Residual 576 653.4 3.787 2.030 308.6
Total 639
P-values were obtained by 4999 permutations of appropriate units, as shown in Table 3 for each term in the
analysis.
P-values in italics were obtained using 4999 Monte Carlo samples from the asymptotic permutation distribution.
Mean squares used in the denominator for each test are also shown in Table 3.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 207
Fig. 5. Means (F1 S.E.) for (a) the total number of fish, (b) the total number of species, (c) the number of P.
scaber (leather jacket) and (d) the number of C. dispilus (demoiselle) in each combination of
YearÂLocationÂHabitat (n = 4 sites per combination of levels and a site consisted of counts summed across
10 transects).
(Sandager’s wrasse), Chironemus marmoratus (hiwihiwi), Parma alboscapularis
(black angelfish), Scorpis violaceus (blue maomao) and Kyphosus sydneyanus (silver
drummer) (Appendix A). All of these fish had strong correlations with canonical
axes, as did Pempheris adspersus (big eye) (Table 7). Average numbers of N.
celidotus also tended to be greater in barrens habitats (Fig. 6a), although this trend
was only statistically significant at Hahei ( P < 0.05), but not at the other locations
(Table 8a).
Some of the common species did not show any strong consistent patterns or differences
across the two habitats, including U. lineatus (goatfish, Fig. 6b, Table 8b), S. lineolatus
(sweep, Fig. 6c, Table 8c) and C. spectabilis (red moki), which was significantly more
208 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 7
Correlations of individual species with the canonical axis for habitat for each of the four locations, as shown in
Fig. 4 (right-hand side)
Name Berghan Home Leigh Hahei Average
Positive correlation (kelp)
Bodianus unimaculatus 0.489 0.639 – 0.232 0.453
Seriola lalandi 0.499 – 0.218 0.561 0.426
Parika scaber 0.681 0.856 0.384 À0.226 0.424
Zeus faber 0.356 0.425 0.475 0.417 0.418
Nemadactylus douglasii 0.491 0.776 0.333 À0.245 0.338
Trachurus novaezelandiae 0.420 0.305 0.250 0.060 0.259
Pseudolabrus miles – 0.270 0.075 0.398 0.247
Chromis dispilus 0.435 0.286 À0.269 0.237 0.172
Odax pullus 0.104 0.456 0.248 À0.172 0.159
Negative correlation (barrens)
Parma alboscapularis À0.409 À0.495 – À0.539 À0.481
Chironemus marmoratus À0.552 À0.023 À0.681 À0.591 À0.462
Coris sandageri À0.352 À0.303 À0.438 À0.551 À0.411
Notolabrus celidotus À0.340 À0.083 À0.380 À0.694 À0.374
Girella tricuspidata À0.414 À0.499 À0.056 À0.442 À0.353
Pempheris adspersus À0.124 À0.102 À0.569 À0.261 À0.264
Notolabrus fucicola À0.119 0.062 À0.028 À0.822 À0.227
Kyphosus sydneyanus À0.125 À0.358 0.218 À0.488 À0.188
Scorpis violaceus 0.347 À0.415 À0.090 À0.456 À0.153
A positive correlation indicates species associated with kelp habitat, while a negative correlation indicates species
associated with barrens habitat. Species are given in decreasing order of the absolute value of their average for the
correlation across the four locations. Dashed lines indicate a species did not occur at that location. Species that
occurred in fewer than 7 transects (out of a total of 640) were not included.
abundant in barrens habitats but only in year 2 and only at Hahei (Fig. 6d, Table 8d, pair-
wise comparisons, P < 0.05).
3.4. Effects of locations
Although locations were originally chosen randomly and so treated as a random
factor in the analyses, it was nevertheless of interest to compare the 4 locations in
terms of the fish assemblages found there, for biogeographic reasons and to compare
results with those of previous studies. Due to its interaction with habitat (Table 3), the
potential differences among locations were considered separately for each habitat. For
each habitat, unconstrained ordinations did not separate assemblages from different
locations very clearly (non-metric MDS plots, Fig. 7, right-hand side), while con-
strained (CAP) ordinations appeared to successfully separate Leigh, Hahei and Home
Point, with assemblages from Berghan Point being less distinct (Fig. 7, right-hand
side). The separation in multivariate space of assemblages from different locations was
slightly more successful for barrens habitats than for kelp forests (Table 9). For barrens
habitats, fish assemblages from each location differed significantly from all other
locations, except for Berghan Point and Home Point, which did not differ significantly
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 209
Fig. 6. Means (F1 S.E.) for (a) the number of N. celidotus (spotty), (b) the number of U. lineatus (goatfish),
(c) the number of S. lineolatus (sweep) and (d) the number of C. spectabilis (red moki) in each combination of
YearÂLocationÂHabitat (n = 4 sites per combination of levels and a site consisted of counts summed across 10
transects).
(Table 10). For kelp forest habitats, fish assemblages did not differ between Berghan
Point and Hahei, but all other comparisons among locations were statistically
significant (Table 10).
The average total number of fish observed per transect was significantly greater at the
two northern locations, Berghan Point and Home Point, than at either Leigh or Hahei (pair-
wise comparisons, P<0.05, Table 6a, Fig. 5a). The greatest average number of species was
found at Hahei for barrens habitats and at Home Point for kelp forests (Table 6b, Fig. 5b).
In year 1, Home Point had a significantly greater average number of species than Leigh or
Hahei, for kelp or barrens habitats. In year 2, this was only true for kelp habitats: in barrens
210 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 8
Results of univariate ANOVAs on selected species
Source df (a) Notolabrus (b) Upeneichthys (c) Scorpis (d) Cheilodactylus
celidotus lineatus lineolatus spectabilis
MS F P MS F P MS F P MS F P
Ye 1 0.046 0.003 0.9584 17.150 0.411 0.5706 228.1 0.398 0.5770 0.693 0.277 0.6396
Lo 3 28.859 6.164 0.0114 97.709 5.018 0.0038 249.9 0.675 0.6062 6.076 3.401 0.0506
Si(Lo) 12 4.682 2.122 0.0126 19.471 2.397 0.0044 370.4 2.275 0.0072 1.787 2.248 0.0104
Ha 1 71.939 3.579 0.1430 15.242 2.524 0.2014 152.0 0.224 0.6660 2.132 0.358 0.5912
YeÂLo 3 15.374 5.017 0.0204 41.756 1.506 0.2226 573.7 1.733 0.2080 2.499 2.860 0.0786
YeÂSi 12 3.064 1.388 0.1678 27.726 3.414 0.0002 331.1 2.034 0.0142 0.874 1.099 0.3466
(Lo)
YeÂHa 1 0.858 0.105 0.8190 4.350 2.096 0.2594 232.9 0.684 0.4786 2.634 0.753 0.4688
LoÂHa 3 20.103 3.849 0.0334 6.040 0.963 0.4586 678.0 1.391 0.2984 5.949 5.654 0.0126
Si(Lo)Â 12 5.223 2.366 0.0064 6.270 0.772 0.7356 487.3 2.994 0.0008 1.052 1.324 0.2068
Ha
YeÂLo 3 8.201 2.050 0.1578 2.076 0.402 0.7866 340.6 0.984 0.4448 3.497 6.829 0.0048
Ha
YeÂSi 12 4.001 1.813 0.0436 5.166 0.636 0.8532 346.3 2.127 0.0068 0.512 0.644 0.8154
(Lo)Ha
Residual 576 2.207 8.123 162.8 0.795
Total 639
P-values were obtained by 4999 permutations of appropriate units, as shown in Table 3 for each term in the
analysis.
P-values in italics were obtained using 4999 Monte Carlo samples from the asymptotic permutation distribution.
Mean squares used in the denominator for each test are also shown in Table 3.
habitats, Hahei had a significantly greater average number of fish species than either
Berghan Point or Leigh (Fig. 5b, pair-wise comparisons).
Significantly greater numbers of P. scaber were found at Hahei than at other
locations for barrens habitats, while for kelp forests, there were significantly greater
average numbers at Home Point compared to the other locations (Fig. 5c, Table 6c).
There were no significant differences among locations in numbers of S. lineolatus, but
C. dispilus was significantly more abundant at the two northern locations: Berghan
Point and Home Point, compared to Leigh and Hahei (Fig. 5d, Table 6d, pair-wise
comparisons, P<0.05). N. celidotus was significantly more abundant at Hahei than at
the other locations in barrens habitats, while there were no significant differences in its
average abundances at different locations for kelp habitats (Table 8a, Fig. 6a, pair-wise
comparisons). U. lineatus was significantly more abundant at Leigh than at the other
locations (Table 8b, Fig. 6b). The only statistically significant differences among
locations for C. spectabilis occurred in year 2 and in barrens habitat only, where
Hahei had significantly greater numbers, on average, than either Leigh or Berghan
Point (Fig. 6d, Table 8d).
Some species were more frequently observed at Home Point than at any other locations,
including P. scaber, P. alboscapularis and N. douglasii (Appendix A). Others had greater
frequencies of occurrence at the two northern locations (Berghan Point and Home Point)
than at the southern locations (Leigh and Hahei), such as C. dispilus, C. sandageri and B.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 211
Fig. 7. Unconstrained metric MDS plots (left) and constrained CAP plots (right) done separately for each habitat
(rows) on the basis of the binomial deviance dissimilarity measure, in each case comparing fish assemblages
among four different locations: Berghan Point, Home Point, Leigh and Hahei. There are eight points for each
combination of the factors, which correspond to the four sites in each of 2 years.
unimaculatus (Appendix A). Alternatively, some species were observed significantly more
frequently at Home Point and Hahei than at any other location: C. spectabilis, G.
tricuspidata, C. marmoratus and Aplodactylus arctidens (marblefish) (Appendix A).
Table 9
Results of CAP analyses examining effects of location within each habitat for 46 species of fish on the basis of the
binomial deviance dissimilarity measure
Habitat m %Var Allocation success (%) d2
B Ho L Ha Total
Barrens 8 91.91 88 63 75 88 78 0.734
Kelp 8 94.56 63 88 50 63 66 0.697
Headings are as defined for Table 5, with B=Berghan Point, Ho=Home Point, L=Leigh, Ha=Hahei.
212 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 10
P-values for NPMANOVA pair-wise comparisons among locations for each Habitat, using 4999 permutations
Habitat Pair-wise P
comparison
Barrens B vs. Ho 0.0618
B vs. L 0.0228
B vs. Ha 0.0106
Ho vs. L 0.0002
Ho vs. Ha 0.0014
L vs. Ha 0.0008
Kelp B vs. Ho 0.0228
B vs. L 0.0050
B vs. Ha 0.0608
Ho vs. L 0.0002
Ho vs. Ha 0.0004
L vs. Ha 0.0006
In each case, the n=10 transects within a site were permuted together as a unit. For each test, there were eight such
units per group (four sites observed per LoÂHa combination in each of 2 years). No corrections have been made
for multiple tests. B=Berghan Point, Ho=Home Point, L=Leigh, Ha=Hahei.
Finally, U. lineatus occurred more frequently at Leigh than anywhere else, and N. fucicola
and T. novaezelandiae occurred more frequently at Hahei than anywhere else (Appendix A).
4. Discussion
There were significant differences in the multivariate structure of fish assemblages
between kelp forest and barrens habitats. These were not independent of species
identity. Differences between habitats included greater frequencies and/or abundances
of P. scaber, C. dispilus, T. novaezelandiae, N. douglasii, B. unimaculatus, O. pullus
and P. miles in kelp forest habitat, while N. celidotus, N. fucicola, G. tricuspidata, C.
sandageri, C. marmoratus, P. alboscapularis, S. violaceus and K. sydneyanus were
more frequent or abundant in barrens habitat. Some of the more common species,
including U. lineatus, S. lineolatus and C. spectabilis showed no strong consistent
effects of these two differing habitats on their distributions. These results differ
somewhat from the results of Choat and Ayling (1987), who focused instead on the
size classes and feeding groups, rather than on individual species, in their description.
The present work examines species’ overall patterns with a broad brush, while the
work by Choat and Ayling (1987) recognized the importance of ontogenetic changes
in diet and habitat use by individual species (e.g., Clements and Choat, 1993; Moran
and Clements, 2002).
One of the most likely reasons for differences in fish assemblages between habitats
is due to depth. Kelp forest and barrens habitats in northeastern New Zealand differ in
their depth distributions (Fig. 2 and see Schiel, 1990). Previous studies have shown
how the numbers and diversity of fish are affected by depth in temperate New
Zealand (e.g., Meekan and Choat, 1997; Brook, 2002) and in tropical coral reefs (e.g.,
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 213
Roberts and Ormond, 1987; Friedlander and Parrish, 1998). We found that depth had
a significant effect on fish assemblage composition (Table 2), although richness and
the total numbers of fish were not significantly related to depth (Table 1). It was clear
in the present study, however, that effects of habitat were not limited to effects of
depth alone (Table 2), suggesting that other aspects of the habitat were important in
structuring assemblages.
It may well be that fundamental biological differences in diet and patterns of foraging
are driving differences in relative frequencies or abundances of species in different
habitats, as suggested by Choat and Ayling (1987). For some species (e.g., O. pullus),
the kelp E. radiata provides a source of food, although Clements and Choat (1993)
suggested that the brown alga Carpophyllum spp. may be preferred by this species. In
contrast, although K. sydneyanus is a herbivore that consumes Ecklonia (Moran and
Clements, 2002), it was found more frequently in barrens habitat in the present study
(Appendix A).
There are many other possible mechanisms that might also explain differences in fish
assemblages in different habitats. There are natural differences in habitat complexity
between kelp and barrens habitats, with kelp forests providing a more complex three-
dimensional biogenic structure. As such, the kelp may provide a refuge from predation,
particularly for juvenile stages (e.g., P. scaber). Also, within kelp forests, we found that
diversity and total abundance of fish were each positively related to variation in kelp
density (Table 1), indicating that the structure and heterogeneity of the habitat may play
a role. Habitat complexity can modify the effects of predation on fish by providing
refuges (e.g., Connell and Jones, 1991; Hixon and Beets, 1993; Caley and St. John,
1996; Beukers and Jones, 1997; Tupper and Boutilier, 1997; Steele, 1999). Alternatively,
differences in kelp density can result in differences in understorey algae and invertebrate
species, which fish may respond to positively or negatively (e.g., Carr, 1989). Connell
(2002) and Wellenreuther and Connell (2002) have demonstrated such effects of habitat-
driven prey availability on spatial patterns of a temperate reef fish (Cheilodactylus
nigripes) from local through to broad scales (hundreds of kilometers). Habitat may also
affect the frequency or intensity of inter-or intra-specific behavioural interactions (e.g.,
Levin et al., 2000).
It is important not to infer too much concerning processes that may be underlying
observed patterns in the results reported here. What we have given is a ‘‘snapshot’’ of
fish we happened to see when and where we saw them. First of all, not all potential
habitats for these fish have been included in this study (e.g., sponge gardens, shallow
areas with dense algal stands of Carpophyllum spp. or other algae, etc.). Thus, for
example, observing increased frequencies of some species in Ecklonia forests should not
be interpreted to mean that this is their ‘‘preferred’’ habitat. Furthermore, we expect
individual species will occur in different habitats at different stages in their life history
or during different behavioural stages (e.g., nesting, feeding, etc.), indeed even during
different stages of the tide or time of day. Quantitative observations we have given here
will have some biological basis, but this will require further and more detailed species-
specific studies.
One important result obtained here is that effects of habitat, although relatively
consistent for the 2 years of observation, did not occur at all locations. In fact,
214 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
multivariate effects of habitat were marginal at Berghan Point ( P = 0.0564) and were not
statistically significant ( P > 0.10) for Leigh (Table 5). This could be one reason that
Choat and Ayling (1987) concluded that habitat effects were independent of species
identity, as many of their survey sites were around the Leigh area. We did find, however,
that where effects of habitat were evident, the direction (i.e., the nature) of multivariate
effects was similar (e.g., Fig. 3). This allowed us to make some generalisations about the
nature of habitat differences, in terms of compositional changes, where they did occur.
We must be cautious, however, in invoking any generalisations of processes that might
have produced observed patterns. For example, Fowler-Walker and Connell (2002) have
suggested that top-down processes are not as important in shaping shallow subtidal
benthic coastal assemblages in South Australia as in New South Wales, because South
Australia lacks the predominance of herbivorous grazers found in regions of New South
Wales. Similarly, variations in effects of habitat among locations within New Zealand
could be due to variation in the abundance or activity of invertebrate grazers and the
relative strength and interactions of bottom-up vs. top-down processes (e.g., Menge,
1992).
Differences in fish assemblages among locations provided interesting biogeographic
information that supports results from previous studies and provided some new insights.
First, there were, on average, greater numbers of species at Home Point than at any
other location, followed by Berghan Point and Hahei, with Leigh, on average, having
the most depauperate fish communities. This is consistent with previous studies of North
Island New Zealand fish fauna, where higher diversity at the Poor Knights Islands, the
Karikari peninsula and Cape Brett was hypothesized to be a consequence of the
influences of the East Auckland Current (e.g., Denham et al., 1984; Choat and Ayling,
1987; Brook, 2002), which does not appear to have much of an influence further south
(i.e., at Leigh). However, in this study we found that, for kelp forest habitats, there was
no significant difference between fish assemblages at Hahei (the southern-most location)
and those at Berghan Point (the northern-most location). So, the influences of the East
Auckland Current appear to extend, at least potentially, to areas as far south as Hahei,
where species such as B. unimaculatus (pigfish), Suezichthys aylingi (crimson cleaner-
fish) and Suezichthys arquatus (rainbowfish) were all recorded in the present study. In
contrast, Leigh may exist as a kind of oceanographic ‘‘backwater’’ of biodiversity for
fish, potentially because of increases in turbidity and/or decreases in relative exposure as
one moves southward into the Hauraki Gulf (Grace, 1983). It is also possible that there
is greater fishing pressure around Leigh, which is not as isolated from human
populations as the other locations in the study. Roberts et al. (1992) also found a
non-linear gradient of fish assemblage structure with latitude in the Red Sea that was
apparently caused by (or at least related to) a negative association of diversity with
turbidity.
The analysis of multivariate variability in fish assemblages at different spatial scales
revealed that the greatest variation occurred at the smallest spatial scale, between
individual transects (Table 4). This is not terribly surprising, given that the spatial
scale of individual transects is not large compared to the high mobility of many fish
species included in these surveys. This result concurs with many studies of inverte-
brates and algae in intertidal and subtidal environments, which have also often found
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 215
the greatest variability to occur at small spatial scales (e.g., Archambault and Bourget,
1996; Underwood and Chapman, 1996; Menconi et al., 1999; Fowler-Walker and
Connell, 2002). In addition, we found that some of the most abundant and ubiquitous
species in the surveys had significant variability from year to year at different sites
and in different habitats (i.e., a significant YearÂSite(Locations)ÂHabitat term for each
of C. dispilus, N. celidotus and S. lineolatus), meaning that predictability at small
scales could prove very difficult for these species. Also, the spatial distribution of
planktivores, such as C. dispilus and S. lineolatus, will likely be heavily dependent on
temporal changes in currents, seasons and tides. We sampled only over a 1-month
period in each of 2 years, but it is clear that variation at different temporal scales
(e.g., tidal, within-day, daily, monthly, etc.) can also be quite large and can pose
difficulties in the assessment of effects of habitat for fish (e.g., Labridae, Connell and
Kingsford, 1998).
Variation from site to site in multivariate assemblages was about the same size as
variability from location to location (Table 4), which concurs with the results of Choat
and Ayling (1987). Observed spatial patterns in fish assemblages were consistent in the
2 years. These results are consistent with previous studies. Jones (1984b) indicated that,
generally, variation in the densities of N. celidotus in different sites was greater than
year-to-year variation over 4 years. Doherty (1987) indicated that, although temporal
variation was detected at virtually every spatial scale found, recruitment was relatively
consistent and hence predictable at geographic and regional scales and Fowler et al.
(1992) found consistent spatial patterns of recruitment from year to year in butterflyfish
on the Great Barrier Reef. In contrast, Tolimieri (1995) found no consistent patterns in
recruitment of damselfish from year to year in the Caribbean and Sale et al. (1984)
found significant year-by-reef interactions in recruitment for nine species of fish on the
Great Barrier Reef. The number of years examined to date (only two) is still too small to
make any inferences about inter-annual variation. Observations of New Zealand’s fish
fauna are needed over longer periods of time (i.e., over several years) and at various
temporal scales to begin to understand potential spatio-temporal interactions or consis-
tencies of patterns.
The observations of patterns given here and by Choat and Ayling (1987) provide a
starting point for further investigations of effects of habitat on distributions of
temperate New Zealand fishes (e.g., Underwood et al., 2000). Experiments, including
manipulations of habitat, observations of feeding, behaviour, competition and preda-
tion, are needed for individual species and groups to further investigate the causal
mechanisms behind any observed relationships or patterns. For example, Jones
(1984a), Carr (1989), Levin and Hay (1996) and Wellenreuther and Connell (2002)
have experimentally altered the presence and/or density of algae to investigate effects
on temperate reef fish. In New Zealand, Jones (1984a) showed experimentally how
juvenile spotties, N. celidotus, recruited in greater abundance in kelp forests and were
negatively affected by decreases in kelp density. Although variation in recruitment
among sites was apparently reflected in adult densities, Jones (1984b) reported
nevertheless, on average, between 4.4 and 7.2 spotties (per 500 m2) in kelp forest
habitat sites, compared to 13.4 and 46.8 (per 500 m2) in shallow broken rock habitat
near Leigh (i.e., Table II therein). Thus, distributions of adult N. celidotus as reported in
216 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Jones (1984b) support the results found here, that is, of greater numbers, on average, in
barrens habitats.
Work along the coast of southern California and reviews of studies in other parts of the
world have emphasised the idea that temperate reef fish, particularly through their
predatory impacts on herbivores, such as urchins, can affect the distributions of kelp
forests (e.g., Cowen, 1983; Holbrook et al., 1990; Tegner and Dayton, 2000). In contrast,
previous studies in New Zealand (including the present study) have emphasised the idea
that temperate reef fishes can respond to changes in habitat, and have not found much
evidence to support the alternative idea that, instead, fish are affecting the structure and
demography of the habitat (e.g., Andrew and Choat, 1982; Choat and Ayling, 1987; Jones,
1988; Schiel, 1990). However, recent studies of the large-scale effects of marine reserves
have revealed that changes in densities of predatory fishes may have strong top-down
effects on community structure (Babcock et al., 1999; Shears and Babcock, 2002). In
particular, increases in predators, such as snapper (Pagrus auratus) and blue cod (Para-
apercis colias), inside of marine reserves resulted in increased rates of predation and
decreases in densities of urchins (E. chloroticus), which in turn has been correlated with
notable increases in the physical extent of kelp forest as opposed to barrens habitat within
reserves (Babcock et al., 1999; Shears and Babcock, 2002). Thus, perhaps not surprisingly,
it would appear that both processes occur: fish populations can affect habitat structure and
habitat, in turn, can affect fish. This recent work, combined with the present study and the
work of Choat and Ayling (1987), suggests that there may be dynamic feedback loops
between indirect effects of fish on habitat structure (especially through predation on
urchins) and direct and indirect effects of habitat structure (particularly the presence and
varying density of kelp forest) on fish communities. However, the relative strength of these
processes may be location-specific, and we must heed warnings of the potential dangers
inherent in generalisations, particularly over large spatial scales (e.g., Underwood and
Petraitis, 1993; Fowler-Walker and Connell, 2002). Clearly, more study is needed and of
more components of the fish fauna in New Zealand to examine potential mechanisms and
the relative importance of these processes.
Acknowledgements
This study could not have been done without the generous help of many diving
assistants and volunteers, including: K. Bloxham, R. Dixon, B. Doak, D. Egli, D. Feary, C.
Galpin, T. Langlois, G. Nesbit, D. Parsons, A. Rapson, M. Rixon, T. Ross-Watt, J.
Saunders, M. Schlegel, N. Tolimieri and T. Willis. We also thank Arthur Cozens, Russ
Babcock, the Leigh Marine Laboratory, Peter Keen and (what was then) the School of
Environmental and Marine Science for logistic support. Some of the data from this study
form part of an investigation into the large-scale patterns of variation in temperate reef
fishes, in collaboration with S. Connell and B. Gillanders (University of Adelaide,
Australian Research Council Discovery Project), who deserve special thanks for all of
their input into the initial experimental design. We also thank K. Clements, S. Connell and
T. Willis for valuable discussions and comments on the manuscript, and N. Shears for
assistance with drawing maps. [AU]
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 217
Appendix A
List of the 46 fish species recorded in the study and their frequencies of occurrence (a)
at each location (out of a possible 160 transects) and (b) in each habitat (out of a possible
320 transects).
Family Species Common (a) Location (b) Habitat
name Berghan Home Leigh Hahei Barrens Kelp
Aplodactylidae Aplodactylus Marblefish 1 12 1 8 14 8
arctidens
Arripidae Arripis trutta Kahawai 0 1 0 5 3 3
Carangidae Decapterus Koheru 2 5 0 2 4 5
koheru
Carangidae Pseudocaranx Trevally 2 6 1 0 6 3
dentex
Carangidae Seriola Kingfish 2 0 1 5 0 8
lalandi
Carangidae Trachurus Jack 15 3 26 42 24 62
novaezelandiae mackerel
Chironemidae Chironemus Hiwihiwi 13 22 7 20 47 15
marmoratus
Dasyatidae Dasyatis Short-tailed 1 5 3 0 3 6
brevicaudata stingray
Dasyatidae Dasyatis Long-tailed 0 3 0 0 1 2
thetidis stingray
Diodontidae Allomycterus Porcupinefish 1 2 0 0 2 1
jaculiferus
Girellidae Girella Parore 13 33 9 21 61 15
tricuspidata
Kyphosidae Kyphosus Silver 2 4 1 4 10 1
sydneyanus drummer
Labridae Anampses Elegant 1 0 0 0 1 0
elegans wrasse
Labridae Bodianus Red pigfish 8 12 0 1 4 17
unimaculatus
Labridae Coris picta Combfish 0 2 0 0 2 0
Labridae Coris Sandager’s 26 27 3 12 44 24
sandageri wrasse
Labridae Notolabrus Spotty 71 69 84 99 175 148
celidotus
Labridae Notolabrus Banded 18 24 4 52 69 29
fucicola wrasse
Labridae Pseudolabrus Orange 1 4 2 0 3 4
luculentus wrasse
Labridae Pseudolabrus Scarlet 0 5 4 9 4 14
miles wrasse
Labridae Suezichthys Rainbowfish 1 0 0 1 1 1
arquatus
(continued on next page)
218 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Appendix A (continued)
Family Species Common (a) Location (b) Habitat
name
Berghan Home Leigh Hahei Barrens Kelp
Labridae Suezichthys Crimson 0 1 0 2 2 1
aylingi cleanerfish
Latridae Cheilodactylus Red moki 62 84 45 70 135 126
spectabilis
Latridae Latridopsis Blue moki 1 0 0 0 0 1
ciliaris
Latridae Latridopsis Copper 0 1 0 1 1 1
forsteri moki
Latridae Nemadactylus Porae 5 11 5 1 2 20
douglasii
Latridae Nemadactylus Tarakihi 0 5 0 0 0 5
macropterus
Monacanthidae Parika scaber Leatherjacket 47 83 28 59 82 135
Mugilidae Aldrichetta Yellow-eyed 0 3 1 0 1 3
forsteri mullet
Mullidae Upeneichthys Goatfish 40 41 70 26 88 89
lineatus
Muraenidae Gymnothorax Yellow moray 0 1 1 2 2 2
prasinus
Myliobatidae Myliobatis Eagle ray 2 6 2 1 5 6
tenuicaudatus
Odacidae Odax pullus Butterfish 4 7 1 8 7 13
Pempheridae Pempheris Big eye 14 9 10 4 21 16
adspersus
Pinguipedidae Parapercis Blue cod 3 5 1 2 8 3
colias
Pomacentridae Chromis Demoiselle 101 115 16 69 137 164
dispilus
Pomacentridae Chromis Yellow 1 0 0 0 0 1
fumea demoiselle
Pomacentridae Chromis Single-spot 1 0 0 1 1 1
hypsilepis demoiselle
Pomacentridae Parma Black 7 21 0 7 35 0
alboscapularis angelfish
Scorpaenidae Scorpaena Northern 0 0 0 2 2 0
cardinalis Scorpionfish
Scorpidae Scorpis Sweep 42 43 48 56 93 96
lineolatus
Scorpidae Scorpis Blue 3 18 2 3 20 6
violaceus maomao
Serranidae Caesioperca Butterfly 0 0 3 0 0 3
lepidoptera perch
Sparidae Pagrus Snapper 19 14 48 3 42 42
auratus
Trachichthyidae Optivus Slender 5 1 2 3 3 8
elongatus roughy
Zeidae Zeus faber John dory 2 1 2 2 0 7
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 219
References
Anderson, M.J., 2001. A new method for non-parametric multivariate analysis of variance. Aust. Ecol. 26,
32 – 46.
Anderson, M.J., 2001. Permutation tests for univariate or multivariate analysis of variance and regression. Can. J.
Fish. Aquat. Sci. 58, 626 – 639.
Anderson, M.J., Robinson, J., 2003. Generalised discriminant analysis based on distances. Aust. N.Z. J. Stat. 45,
301 – 318.
Anderson, M.J., ter Braak, C.J.F., 2003. Permutation tests for multi-factorial analysis of variance. J. Stat. Comput.
Simul. 73, 85 – 113.
Anderson, M.J., Willis, T.J., 2003. Canonical analysis of principal coordinates: a useful method of constrained
ordination for ecology. Ecology 84, 511 – 525.
Andrew, N.L., Choat, J.H., 1982. The influence of predation and conspecific adults on the abundance of juvenile
Evechinus chloroticus (Echinoidea: Echinodermata). Oecologia 54, 80 – 87.
Andrew, N.L., Mapstone, B.D., 1987. Sampling and the description of spatial pattern in marine ecology. Ocean-
ogr. Mar. Biol., Ann. Rev. 25, 39 – 90.
Archambault, P., Bourget, E., 1996. Scales of coastal heterogeneity and benthic intertidal species richness,
diversity and abundance. Mar. Ecol. Prog. Ser. 136, 111 – 121.
Ault, T.R., Johnson, C.R., 1998. Spatially and temporally predictable fish communities on coral reefs. Ecol.
Monogr. 68, 25 – 50.
Babcock, R.C., Kelly, S., Shears, N.T., Walker, J.W., Willis, T.J., 1999. Changes in community structure in
temperate marine reserves. Mar. Ecol. Prog. Ser. 189, 125 – 134.
Benedetti-Cecchi, L., 2001. Variability in abundance of algae and invertebrates at different spatial scales on rocky
sea shores. Mar. Ecol. Prog. Ser. 215, 79 – 92.
Beukers, J.S., Jones, G.P., 1997. Habitat complexity modifies the impact of piscivores on a coral reef fish
population. Oecologia 114, 50 – 59.
Brock, R.E., 1982. A critique of the visual census method for assessing coral reef fish populations. Bull. Mar. Sci.
32, 269 – 276.
Brook, F.J., 2002. Biogeography of near-shore reef fishes in northern New Zealand. J. R. Soc. N.Z. 32, 243 – 274.
Caley, M.J., St. John, J., 1996. Refuge availability structures assemblages of tropical reef fishes. J. Anim. Ecol.
65, 414 – 428.
Cao, Y., Williams, W.P., Bark, A.W., 1997. Similarity measure bias in river benthic Aufwuchs community
analysis. Water Environ. Res. 69 (1), 95 – 106.
Carr, M.H., 1989. Effects of macroalgal assemblages on the recruitment of temperate zone reef fishes. J. Exp.
Mar. Biol. Ecol. 126, 59 – 76.
Chesson, P., 1998. Spatial scale in the study of reef fishes: a theoretical perspective. Aust. J. Ecol. 23,
209 – 215.
Choat, J.H., Ayling, A.M., 1987. The relationship between habitat structure and fish faunas on New Zealand
reefs. J. Exp. Mar. Biol. Ecol. 110, 257 – 284.
Clark, P.J., 1952. An extension of the coefficient of divergence for use with multiple characters. Copeia 1952,
61 – 64.
Clarke, K.R., 1993. Non-parametric multivariate analyses of changes in community structure. Aust. J. Ecol. 18,
117 – 143.
Clements, K.D., Choat, J.H., 1993. Influence of season, ontogeny and tide on the diet of the temperate marine
herbivorous fish Odax pullus (Odacidae). Mar. Biol. 117, 213 – 220.
Connell, S.D., 2002. Effects of predator and prey on a foraging reef fish: implications for understanding density-
dependent growth. J. Fish Biol. 60, 1551 – 1561.
Connell, S.D., Jones, G.P., 1991. The influence of habitat complexity on postrecruitment processes in a temperate
reef fish population. J. Exp. Mar. Biol. Ecol. 151, 271 – 294.
Connell, S.D., Kingsford, M.J., 1998. Spatial, temporal and habitat-related variation in the abundance of large
predatory fish at One Tree Reef, Australia. Coral Reefs 17, 49 – 57.
Cowen, R.K., 1983. The effect of sheephead (Semicossyphus pulcher) predation on red sea urchin (Strong-
gylocentrotus franciscanus) populations: an experimental analysis. Oecologia 58, 249 – 255.
220 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
DeMartini, E.E., Roberts, D.A., 1990. The effects of giant kelp (Macrocystis) on the density and abundance of
fishes in a cobble-bottom kelp forest. Bull. Mar. Sci. 46, 287 – 300.
Denham, R.N., Bannister, K.M., Guthrie, K.M., Crook, F.G., 1984. Surveys of the east Auckland and east cape
currents, New Zealand. Aust. J. Mar. Freshw. Res. 35, 491 – 504.
Doherty, P.J., 1987. The replenishment of populations of coral reef fishes, recruitment surveys, and the problems
of variability manifest on multiple scales. Bull. Mar. Sci. 41, 411 – 422.
Fowler, A.J., Doherty, P.J., Williams, D.MCB., 1992. Multi-scale analysis of recruitment of a coral reef fish on
the Great Barrier Reef. Mar. Ecol. Prog. Ser. 82, 131 – 141.
Fowler-Walker, M.J., Connell, S.D., 2002. Opposing states of subtidal habitat across temperate Australia: con-
sistency and predictability in kelp canopy-benthic associations. Mar. Ecol. Prog. Ser. 240, 49 – 56.
Friedlander, A.M., Parrish, J.D., 1998. Habitat characteristics affecting fish assemblages on a Hawaiian coral reef.
J. Exp. Mar. Biol. Ecol. 224, 1 – 30.
´ ´
Garcıa-Charton, J.A., Perez-Ruzafa, A., 2001. Spatial pattern and the habitat structure of a Mediterranean reef
fish local assemblage. Mar. Biol. 138, 917 – 934.
Gower, J.C., Legendre, P., 1986. Metric and Euclidean properties of dissimilarity coefficients. J. Classif. 3, 5 – 48.
Grace, R.V., 1983. Zonation of sublittoral rocky bottom marine life and its changes from the outer to the inner
Haruaki Gulf, north-eastern New Zealand. Tane 29, 97 – 108.
Hixon, M.A., Beets, J.P., 1993. Predation, prey refuges and the structure of coral-reef fish assemblages. Ecol.
Monogr. 63, 77 – 101.
Holbrook, S.J., Schmitt, R.J., Ambrose, R.F., 1990. Biogenic habitat structure and characteristics of temperate
reef fish assemblages. Aust. J. Ecol. 15, 489 – 503.
Holbrook, S.J., Forrester, G.E., Schmitt, R.J., 2000. Spatial patterns in abundance of a damselfish reflect
availability of suitable habitat. Oecologia 122, 109 – 120.
Hurlbert, S.H., 1984. Pseudoreplication and the design of ecological field experiments. Ecol. Monogr. 54,
187 – 211.
Jones, G.P., 1984. Population ecology of the temperate reef fish Pseudolabrus celidotus Bloch and Schneider
(Pisces: Labridae): I. Factors influencing recruitment. J. Exp. Mar. Biol. Ecol. 75, 257 – 276.
Jones, G.P., 1984. Population ecology of the temperate reef fish Pseudolabrus celidotus Bloch and Schneider
(Pisces: Labridae): II. Factors influencing adult density. J. Exp. Mar. Biol. Ecol. 75, 277 – 303.
Jones, G.P., 1988. Ecology of rocky reef fish of north-eastern New Zealand: a review. N.Z. J. Mar. Freshw. Res.
22, 445 – 462.
Kotliar, N.B., Wiens, J.A., 1990. Multiple scales of patchiness and patch structure: a hierarchical framework for
the study of heterogeneity. Oikos 59, 253 – 260.
Lachenbruch, P.A., Mickey, M.R., 1968. Estimation of error rates in discriminant analysis. Technometrics 10, 1 – 11.
Legendre, P., Fortin, M.J., 1989. Spatial pattern and ecological analysis. Vegetatio 80, 107 – 138.
Legendre, P., Legendre, L., 1998. Numerical Ecology, 2nd English edition Elsevier, Amsterdam.
Levin, P.S., Hay, M.E., 1996. Responses of temperate reef fishes to alterations in algal structure and species
composition. Mar. Ecol. Prog. Ser. 134, 37 – 47.
Levin, P.S., Toimieri, N., Nicklin, M., Sale, P.F., 2000. Integrating individual behavior and population ecology:
the potential for habitat-dependent population regulation in a reef fish. Behav. Ecol. 11, 565 – 571.
Lincoln-Smith, M.P., 1989. Improving multispecies rocky reef fish censuses by counting different groups of
species using different procedures. Environ. Biol. Fishes 26, 29 – 37.
Manly, B.F.J., 1997. Randomization, Bootstrap and Monte Carlo Methods in Biology Chapman & Hall, London.
McArdle, B.H., Anderson, M.J., 2001. Fitting multivariate models to community data: a comment on distance-
based redundancy analysis. Ecology 82, 290 – 297.
McClanahan, T.R., Arthur, R., 2001. The effect of marine reserves and habitat on populations of east African
coral reef fishes. Ecol. Appl. 11, 559 – 569.
Meekan, M.G., Choat, J.H., 1997. Latitudinal variation in abundance of herbivorous fishes: a comparison of
temperate and tropical fishes. Mar. Biol. 128, 373 – 383.
Menconi, M., Benedetti-Cecchi, L., Cinelli, F., 1999. Spatial and temporal variability in the distribution of algae
and invertebrates on rocky shores in the northwest Mediterranean. J. Exp. Mar. Biol. Ecol. 233, 1 – 23.
Menge, B.A., 1992. Community regulation: under what conditions are bottom-up factors important on rocky
shores?. Ecology 73, 755 – 765.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 221
Moran, D., Clements, K.D., 2002. Diet and endogenous carbohydrases in the temperate marine herbivorous fish
Kyphosus sydneyanus. J. Fish Biol. 60, 1190 – 1203.
Paulin, C., Stewart, A., Roberts, C., McMillan, P., 1989. New Zealand fish: a complete guide. National Museum
of New Zealand Miscellaneous Series, vol. 19. GP Books, Wellington, New Zealand, pp. 1 – 279.
Roberts, C.M., Ormond, R.F.G., 1987. Habitat complexity and coral reef fish diversity and abundance on Red Sea
fringing reefs. Mar. Ecol. Prog. Ser. 41, 1 – 8.
Roberts, C.M., Dawson Shepard, A.R., Ormond, R.F.G., 1992. Large-scale variation in assemblage structure of
Red Sea butterflyfishes and angelfishes. J. Biogeogr. 19, 239 – 250.
Sale, P.F., 1998. Appropriate spatial scales for studies of reef-fish ecology. Aust. J. Ecol. 23, 202 – 208.
Sale, P.F., Sharp, B.J., 1983. Corrections for bias in visual transect censuses of coral reef fishes. Coral Reefs 2,
37 – 42.
Sale, P.F., Doherty, P.J., Eckert, G.J., Douglas, W.A., Ferrell, D.J., 1984. Large scale spatial and temporal
variation in recruitment of fish populations to fish populations on coral reefs. Oecologia 64, 191 – 198.
Schiel, D.R., 1990. Macroalgal assemblages in New Zealand: structure, interactions and demography. Hydro-
biologia 192, 59 – 76.
Searle, S.R., Casella, G., McCulloch, C.E., 1992. Variance Components Wiley, Toronto.
Shears, N.T., Babcock, R.C., 2002. Marine reserves demonstrate top-down control of community structure on
temperate reefs. Oecologia 132, 131 – 142.
Snedecor, G.W., 1946. Statistical Methods, 4th edition Iowa State College Press, Ames, Iowa.
Steele, M.A., 1999. Effects of shelter and predators on reef fishes. J. Exp. Mar. Biol. Ecol. 233, 65 – 79.
Stephens Jr., J.S., Morris, P.A., Zerba, K., Love, M., 1984. Factors affecting fish diversity on a temperate reef: the
fish assemblage of Palos Verdes Point, 1974 – 1981. Environ. Biol. Fishes 11, 259 – 275.
Tegner, M.J., Dayton, P.K., 2000. Ecosystem effects of fishing in kelp forest communities. ICES J. Mar. Sci. 57,
579 – 589.
Thompson, A.A., Mapstone, B.D., 1997. Observer effects and training in underwater visual surveys of reef fishes.
Mar. Ecol. Prog. Ser. 154, 53 – 63.
Tolimieri, N., 1995. Effects of microhabitat characteristics on the settlement and recruitment of a coral reef fish at
two spatial scales. Oecologia 102, 52 – 63.
Tolimieri, N., 1998. Constrasting effects of microhabitat use on large-scale adult abundance in two families of
Caribbean reef fishes. Mar. Ecol. Prog. Ser. 167, 227 – 239.
Tolimieri, N., 1998. The relationship among microhabitat characteristics, recruitment and adult abundance in the
stoplight parrotfish, Sparisoma viride, at three spatial scales. Bull. Mar. Sci. 62, 253 – 268.
Tupper, M., Boutilier, R.G., 1997. Effects of habitat on settlement, growth, predation risk and survival of a
temperate reef fish. Mar. Ecol. Prog. Ser. 151, 225 – 236.
Underwood, A.J., Chapman, M.G., 1996. Scales of spatial patterns of distribution of intertidal invertebrates.
Oecologia 107, 212 – 224.
Underwood, A.J., Chapman, M.G., 1998. A method for analyzing spatial scales of variation in composition of
assemblages. Oecologia 117, 570 – 578.
Underwood, A.J., Petraitis, P.S., 1993. Structure of intertidal assemblages in different locations: how can local
processes be compared?. In: Ricklefs, R., Schluter, D. (Eds.), Species Diversity in Ecological Communities.
University of Chicago Press, Chicago, pp. 38 – 51.
Underwood, A.J., Chapman, M.G., Connell, S.D., 2000. Observations in ecology: you can’t make progress on
processes without understanding the patterns. J. Exp. Mar. Biol. Ecol. 250, 97 – 115.
Wellenreuther, M., Connell, S.D., 2002. Response of predators to prey abundance: separating the effects of prey
density and patch size. J. Exp. Mar. Biol. Ecol. 273, 61 – 71.
Wiens, J.A., 1989. Spatial scaling in ecology. Funct. Ecol. 3, 385 – 397.
Willis, T.J., 2001. Visual census methods underestimate density and diversity of cryptic reef fishes. J. Fish Biol.
59, 1408 – 1411.
Willis, T.J., Millar, R.B., Babcock, R.C., 2000. Detection of spatial variability in relative density of fishes:
comparison of visual census, angling, and baited underwater video. Mar. Ecol. Prog. Ser. 198, 249 – 260.
Yates, F., Cochran, W.G., 1938. The analysis of groups of experiments. J. Agric. Sci. 28, 556 – 580.
305 (2004) 191 – 221
www.elsevier.com/locate/jembe
Spatial variation and effects of habitat on temperate
reef fish assemblages in northeastern New Zealand
Marti J. Anderson *, Russell B. Millar
Department of Statistics, Tamaki Campus, University of Auckland, Private Bag 92019, Auckland, New Zealand
Received 15 January 2003; received in revised form 5 November 2003; accepted 22 December 2003
Abstract
Reef-associated fishes can respond to changes in habitat structure and the nature of their
response can change with different spatial scales of observation. A structured hierarchical
mensurative sampling design was used to sample temperate reef fish assemblages in northeastern
New Zealand at several spatial scales over 2 years. The three spatial scales examined were tens
of meters (transects), hundreds to thousands of meters (sites) and hundreds of kilometers
(locations). We tested the hypothesis that fish assemblages differed between kelp forest habitat
(relatively dense stands of the kelp, Ecklonia radiata (C. Agardh) J. Agardh, median depth=13.5
m) and barrens habitat (rocky reef dominated by turfing and encrusting red algae and the
grazing urchin, Evechinus chloroticus (Valenciennes), median depth=6.7 m). Recently developed
multivariate techniques were used to test for and quantify multivariate variation at different
spatial scales. There were significant effects of habitat on the spatial distribution of fish
assemblages, characterised by greater abundances or frequencies of Parika scaber, Chromis
dispilus, Trachurus novaezelandiae, Nemadactylus douglasii, Bodianus unimaculatus, Odax
pullus and Pseudolabrus miles in kelp forest habitat, and greater abundances or frequencies of
Notolabrus celidotus, Notolabrus fucicola, Girella tricuspidata, Coris sandageri, Chironemus
marmoratus, Parma alboscapularis, Scorpis violaceus and Kyphosus sydneyanus in barrens
habitat. Some of the more common species, including Upeneichthys lineatus, Scorpis lineolatus
and Cheilodactylus spectabilis showed no strong consistent effects of these two differing habitats
on their distributions. There was, however, a significant HabitatÂLocations interaction: effects of
habitat did not occur at all locations. Variability was highest at the scale of individual transects
and variability from site to site and from location to location was comparable. Spatial variation
was large compared to inter-annual variation, which was minimal, and spatial patterns were
consistent in the 2 years examined. Further experiments, including manipulations, are required to
* Corresponding author. Tel.: +64-9-373-7599x85052; fax: +64-9-373-7000.
E-mail address: mja@stat.auckland.ac.nz (M.J. Anderson).
0022-0981/$ - see front matter D 2004 Elsevier B.V. All rights reserved.
doi:10.1016/j.jembe.2003.12.011
192 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
understand what mechanisms and processes might be driving these patterns. This study, coupled
with results from previous studies, suggests that there may be a dynamic inter-play between
effects of habitat on fish and effects of fish on biogenic habitat, such as kelp forests.
D 2004 Elsevier B.V. All rights reserved.
Keywords: Hierarchical experimental design; Kelp forest; Multivariate analysis; Temperate reef fish; Urchin
grazing; Variance components
1. Introduction
An important goal in ecology is to understand patterns of distributions of organisms
by reference to the habitat available to them in the environment. For reef-associated
fishes, there is abundant evidence that the structure of the habitat has important effects
on spatial distributions of populations, both in tropical coral reefs (Roberts and
Ormond, 1987; Tolimieri, 1995; Caley and St. John, 1996; Friedlander and Parrish,
1998; Tolimieri, 1998a; Holbrook et al., 2000; McClanahan and Arthur, 2001) and in
temperate rocky reef systems (Choat and Ayling, 1987; Jones, 1988; Holbrook et al.,
1990; Connell and Jones, 1991; Carr, 1989; Levin and Hay, 1996; Tupper and
´ ´
Boutilier, 1997; Garcıa-Charton and Perez-Ruzafa, 2001). For example, for temperate
reefs in northeastern New Zealand, Choat and Ayling (1987) described two distinct
assemblages of fishes associated with two different habitat types: (i) areas dominated
by grazing urchins (Evechinus chloroticus (Valenciennes)), with cover by encrusting
and turfing red algae (called ‘‘barrens’’ habitat) and (ii) areas of relatively dense
stands of the laminarian kelp Ecklonia radiata (C. Agardh) J. Agardh (referred to as
‘‘kelp forests’’). They found that small wrasses were more abundant in the shelter of
kelp forests and increased with increases in kelp density, while larger carnivores were
more abundant in urchin-grazed areas (Choat and Ayling, 1987). Jones (1984a) found
that juveniles of the wrasse, Notolabrus celidotus (spotty), were positively associated
with the density of macroalgae on New Zealand reefs. Furthermore, densities of
recruits decreased with the experimental removal of kelp and increased with the
experimental addition of kelp (through removal of urchins, Jones, 1984a). The
consequences of the use of kelp forests by younger fish for predicting distributions
of adult populations are, however, largely unknown.
There is some evidence for the effects of kelp habitat on reef fishes in other parts
of the world. However, Holbrook et al. (1990) found only weak differences in species
composition among reefs with different types of algal habitats in Southern California,
and mixed results were obtained by other workers in this region (e.g., Stephens et al.,
1984; Carr, 1989; DeMartini and Roberts, 1990). In addition, in New Zealand, Choat
and Ayling (1987) suggested that the differences found between habitats in their study
were independent of species identity, being driven instead by differences in the
biology of feeding preferences of fishes at different ontogenetic stages. Thus, it is
unclear whether the species composition and relative abundances of individual fish
species differ in consistent or predictable ways between these two identifiable habitats
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 193
in northeastern New Zealand. A structured mensurative experiment is needed to
investigate this.
A further important aspect of understanding spatial distributions with respect to
habitat characteristics is to recognize that observed patterns are dependent on the
spatial scale of observation (e.g., Andrew and Mapstone, 1987; Wiens, 1989; Tolimieri,
1995; Chesson, 1998; Sale, 1998). There have been several studies of spatial patterns
of distributions of fish at several spatial scales, from meters up to hundreds or
thousands of kilometers (e.g., Choat and Ayling, 1987; Doherty, 1987; Fowler et al.,
1992; Tolimieri, 1998a; Ault and Johnson, 1998; Connell and Kingsford, 1998; Garcıa- ´
´
Charton and Perez-Ruzafa, 2001; Connell, 2002). Effects of habitat at small spatial
scales can provide predictive power at larger spatial scales for some coral reef fish
species (i.e., damselfish, Tolimieri, 1995; Holbrook et al., 2000), but not for others
(i.e., stoplight parrotfish, Tolimieri, 1998a,b). Jones (1988) suggested that it is the
comparison of the magnitude of variation at each spatial scale of interest that will
provide the necessary framework for formulating hypotheses about relevant processes
for reef fishes. Hierarchical spatially structured sampling programs provide a means of
partitioning and quantifying the magnitude of variation at different spatial scales
(Andrew and Mapstone, 1987; Underwood and Chapman, 1996; Underwood et al.,
2000). Indeed, it may not be possible to understand effects of habitat, which is a
spatial phenomenon, involving patchiness and heterogeneity, without an understanding
of variability at different spatial scales (Kotliar and Wiens, 1990; Underwood and
Chapman, 1996; Underwood et al., 2000). For example, Fowler-Walker and Connell
(2002) demonstrated scale-dependent associations of understorey algae in E. radiata
forests, with weak or variable local-scale responses to habitat, but strong and consistent
patterns at regional scales (thousands of kilometers). Although some explicit measure-
ments of variation in effects of habitat on fish species at different spatial scales have
been made by certain workers (e.g., Tolimieri, 1998a; Doherty, 1987; Fowler et al.,
´ ´
1992; Garcıa-Charton and Perez-Ruzafa, 2001), this has not generally been done using
structured hierarchical designs, such as those used to advantage for benthic invertebrate
and algal assemblages (e.g., Underwood and Chapman, 1996; Menconi et al., 1999;
Benedetti-Cecchi, 2001; Fowler-Walker and Connell, 2002), but see Connell and
Kingsford (1998) and Connell (2002).
One stumbling block towards measuring spatial variation in fish or other kinds of
assemblages has been the lack of an appropriate method for assessing multivariate
variation for several species simultaneously, that is, to measure and quantify variation
in multivariate assemblages in each scale of a hierarchy of spatial scales. Available
distance-based multivariate methods that are realistic for non-normal counts of species
abundances are either unable to obtain independent partitions of multivariate variation for
such complex designs (such as ANOSIM, Clarke, 1993; or Mantel correlograms, Legendre
and Fortin, 1989), or require large numbers of replicates to avoid problems of non-
independence (Underwood and Chapman, 1998).
Recent developments in non-parametric multivariate analysis provide a method for
analyzing multivariate assemblages on the basis of any distance or dissimilarity
measure, while also allowing the data to be partitioned according to any experimental
design, including nested hierarchies (Anderson, 2001a; McArdle and Anderson, 2001).
194 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Non-parametric multivariate analysis of variance (NPMANOVA, Anderson, 2001a) can
be used to partition variability and to provide measures of multivariate variability at
different scales in a structured hierarchical design, in a manner directly analogous to
univariate partitioning using ANOVA (e.g., Searle et al., 1992; Benedetti-Cecchi,
2001). The statistical significance of multivariate variance components can also be
tested using appropriate methods of permutation (Anderson, 2001b; Anderson and ter
Braak, 2003).
The purpose of the present investigation was to conduct a mensurative, observa-
tional experiment (sensu Hurlbert, 1984) to examine the potential effects of habitat on
fish assemblages in northeastern New Zealand at several spatial scales. Before
hypotheses concerning any potential emerging patterns or processes can be developed,
good observational data and quantitative measures of spatial variability in fish
assemblages are extremely useful—for individual species and for the multivariate
assemblage as a whole (Underwood et al., 2000). More specifically, we wished to test
the null hypothesis (following Choat and Ayling, 1987) that there are no differences in
the composition of fish assemblages in kelp forest habitat vs. barrens habitat on
subtidal rocky reefs (i.e., that differences are ‘‘independent of species identity’’). To
determine the generality of any potential patterns, effects of habitat need to be
examined at several locations and at several times (e.g., Yates and Cochran, 1938;
Snedecor, 1946; Underwood and Petraitis, 1993). Furthermore, as effects of habitat on
fish may or may not ‘‘scale up’’ to provide good prediction at larger spatial scales
(e.g., Tolimieri, 1998a), depending on the species (e.g., Ault and Johnson, 1998), we
predicted that the multivariate effects of habitat would interact with variability at
different spatial scales and that different species would show different patterns in this
regard. We also tested the hypothesis that there is a significant relationship between
fish assemblage structure and the density of kelp forests, as a relationship with kelp
density was observed for total numbers of large carnivorous fish by Choat and Ayling
(1987).
The spatial scales of investigation included in the study were tens of meters
(transects), hundreds of meters to kilometers (sites) and hundreds of kilometers
(locations). These large-scale surveys were performed in each of 2 years using a
structured hierarchical experimental design to investigate effects of habitat on reef
fishes. This is the first study, to our knowledge, that uses multivariate methods to obtain
independent quantitative measures of variability in fish assemblages at different spatial
scales.
In northeastern New Zealand, it is known that patches of kelp forest habitat tend to
occur at deeper depths than patches of barrens habitat, although there is some overlap in
their natural depth ranges (e.g., Schiel, 1990). Brook (2002) found a significant positive
relationship between depth (up to 45 m) and species richness of fishes, while depth-
stratified sampling by Meekan and Choat (1997) showed significant differences in
abundances of several prominent herbivorous fishes among different depth strata. It is
not known, however, the extent to which effects of depth may mediate effects of habitat.
The present study therefore also measured depth as a covariable for analyses, to test the
hypothesis that significant differences in assemblages due to habitat, if present, could not
be fully explained by a relationship of assemblages with depth.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 195
2. Methods
2.1. Underwater visual sampling of fish
The study used a structured hierarchical experimental design. Four locations, separated
by hundreds of kilometers, were sampled along the northeastern coast of New Zealand
(Fig. 1). These included, from South to North, Hahei (36j50.86VS, 175j49.32VE), Leigh
(36j17.43VS, 174j48.82VE), Home Point (35j19.38VS, 174j21.38VE) and Berghan Point
(34j55.78VS, 173j32.72VE) (Fig. 1). Within each location, fish were sampled at each of
four different randomly located sites, separated by hundreds to thousands of meters. At
each site, two habitats were investigated: kelp forests (i.e., areas characterised by relatively
dense cover of the kelp E. radiata) and ‘‘barrens’’ (i.e., areas characterised by little or no
macroalgal cover and dominated by the grazing urchin E. chloroticus). The range of
depths was 2 –20 m, although the majority of the observations were in the range from 5 to
15 m.
Within each habitat, divers on SCUBA did a visual survey by swimming along a
transect and identifying and recording the number of each species of fish observed
within a distance of 2.5 m on either side of the transect. Taxonomic authorities for all
fish species named herein are available in Paulin et al. (1989). All fish seen were
recorded except for notoriously cryptic or very small species (e.g., Tripterygiidae and
Gobiidae), which were not included, as we could not be confident that they were being
reliably seen over these spatial scales (e.g., Lincoln-Smith, 1989; Willis, 2001). Each
diver carried a tape measure and swam a ‘‘run-in’’ distance of 5 m before beginning the
actual survey for each transect, which was 25 m long after the run-in. Thus the total area
sampled per transect was 5Â25 m. There were n=10 transects sampled haphazardly
within each habitat at each site. Although there are some known biases inherent in using
the method of underwater visual surveys to count fish (e.g., Brock, 1982; Sale and
Sharp, 1983; Lincoln-Smith, 1989; Thompson and Mapstone, 1997; Willis et al., 2000),
this same sampling methodology was used across the entire experimental design. This
allows valid comparisons to be made across habitats and across different spatial scales,
even though unbiased estimates of population densities may not be obtained for some
fish species from these data. It was not possible logistically for a single diver to perform
the entire sampling design, so several divers participated in the study. There was likely
variation due to different observers, although this was not explicitly measured. This was
not, however, related in any systematic fashion with differences across habitats because
each diver involved in the study sampled both types of habitat. Divers also recorded the
depth at the beginning and at the end of each transect (in meters), which generally
differed by no more than 1 – 2 m. The average depth measured for each transect was
used for analyses. Also, in the first year of sampling for kelp habitats, the density of
kelp in the forests was estimated. This was done by a second diver swimming behind
the diver counting fish. A 1-m2 quadrat was placed alongside the transect tape (used by
the diver counting fish) at each of five positions along the tape (approximately 5, 10, 15,
20 and 25 m) and the number of kelp plants of the species E. radiata per quadrat was
recorded. The average and the standard deviation for kelp density were then calculated
for each transect based on these five observations.
196 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Fig. 1. Map of northeastern New Zealand showing the four locations for the study and the sites where surveys
were done. Grey circles indicate kelp habitat and black triangles indicate barrens habitat sampled at each site.
The sampling was done within a 1-month period in each of 2 consecutive years during
the southern hemisphere’s (austral) summer: from 30 November to 21 December 2000
(year 1) and from 7 January to 5 February 2002 (year 2). The locations were not sampled
consecutively from north to south or vice versa, so as not to confound the latitudinal
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 197
gradient (if any) with time of sampling within the month. In addition, all sites were located
outside of any established marine reserves (i.e., at Hahei and Leigh). GPS coordinates
were recorded at each site in year 1 and the same sites (although not the same transects)
were re-visited in year 2 by reference to these GPS coordinates. Due to the high mobility
of fish, the 2 years were considered as independent samples.
2.2. Multivariate statistical methods
The experimental design therefore consisted of four factors: year (two levels, fixed),
location (four levels, random), site (four levels, random, nested in location) and habitat
(two levels, fixed, crossed with all other factors). Thus, with n=10 transects, there were a
total of 640 observation units in the data set. There were 46 species of fish (variables)
recorded in the study that were included in multivariate analyses (listed in Appendix A).
Although it would have been logical to treat the factor ‘‘year’’ as random in the present
design, with only 2 years of data, inferences concerning inter-annual variation in general
could not be viewed as very precise. Thus, ‘‘year’’ was treated as a fixed factor, which
focused our analyses on results obtained for these 2 years only. With a greater number of
years of sampling anticipated in the future, we would treat this as a random effect.
The fish species variables were highly skewed and contained a great many zeros,
making traditional analyses (which assume normality of errors) unsuitable, so non-
parametric approaches were used. NPMANOVA (Anderson, 2001a; McArdle and Ander-
son, 2001) was used to analyse the multivariate data set in response to the complete
experimental design (including interactions). This method allows multivariate data to be
analysed on the basis of any distance or dissimilarity measure of choice, in response to any
complex experimental design, with P-values obtained using permutations. We used scale-
invariant binomial deviance as a dissimilarity measure (described below) for the NPMA-
NOVA. The tests do not assume that the original variables conform to a multivariate
normal or to any other specified distribution. Although the approach does not explicitly
assume common variances among groups, it will (like ANOSIM, Clarke, 1993) be
sensitive to differences in multivariate dispersion. In addition, the results of NPMANOVA
were also used to estimate the sizes of multivariate pseudo variance components from the
analyses, relying on the analogy of NPMANOVA as an ANOVA based on distances and
univariate ANOVA estimators based on mean squares (e.g., Searle et al., 1992). Doing this
(and indeed applying NPMANOVA to a complex multi-factorial design in general) does
require us to assume that the linear ANOVA model can be applied successfully to achieve
a partitioning of the squared inter-point dissimilarities for inferences to be made about the
multivariate assemblage. We feel this is reasonable provided the choice of dissimilarity
measure accurately reflects relevant qualities of the multivariate assemblage of interest.
For each term in the analysis, 4999 permutations of raw data units were done to obtain
P-values (e.g., Manly, 1997). Care was taken to ensure that the correct permutable units
were used to obtain a valid permutation test of each term in the analysis (Anderson,
2001b; Anderson and ter Braak, 2003). For example, to test the following term:
YearÂLocationÂHabitat, the 64 cell units corresponding to the 64 combinations of
YearÂSite(Location)ÂHabitat were permuted (i.e., the 10 transects within each of these
cells were kept together as a unit). This is because the denominator mean square for the
198 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
test of YearÂLocationÂHabitat is the mean square for YearÂSite(Location)ÂHabitat,
thus it is these units (indicated by the denominator of the F-ratio) that are appropriate to
permute under the null hypothesis (Anderson, 2001b; Anderson and ter Braak, 2003). An
important assumption underlying each test is therefore the exchangeability of the correct
errors (for each particular term being tested in the analysis) under each null hypothesis of
interest. In some cases, this restriction on the permutable units meant that there were not
enough possible permutations to get a reasonable test. For these situations, a Monte Carlo
sample was drawn from the theoretical asymptotic permutation distribution (Anderson
and Robinson, 2003). We re-iterate that these permutation tests assume not just
exchangeability, but also that a linear model on dissimilarities is appropriate for choosing
reasonable test-statistics for testing multivariate hypotheses in an ANOVA framework,
analogous to those used in univariate ANOVA. All analyses were done using specialised
software, written in FORTRAN.
To compare whole fish assemblages to quantitative variables (i.e., depth, average kelp
density and the standard deviation of kelp density), non-parametric multivariate multiple
regression was used on the basis of the binomial deviance dissimilarity measure, using
4999 random permutations (McArdle and Anderson, 2001). We were particularly
concerned to test for any significant effects of habitat over and above the potential effect
of depth, as kelp forests and barrens habitats occur naturally at different depths, on
average. Thus, the effect of habitat was tested, with depth as a covariable in the analysis,
using non-parametric multivariate multiple regression and 4999 permutations of residuals
under the reduced model (e.g., Anderson, 2001b). We also used ordinary least squares
regression to compare each of the three quantitative variables to the total number of
species per transect and the total number of individuals (the latter after transforming to
ln(x+1)). Individual transects were removed from the analysis whenever the quantitative
predictor variable of interest was missing. Also, for these and any other analyses of the
total number of fish (as a univariate variable), the following schooling species that can
occur sporadically in clumps of thousands were not included: Trachurus novaezelandiae
(jack mackerel), Aldrichetta forsteri (yellow-eyed mullet) and Decapterus koheru
(koheru). Far from being dominant species, mullet and koheru, for example, were only
observed in 4 and 16 out of 640 transects, respectively. More importantly, these species
were excluded from sums of the total number of individuals primarily because no certainty
regarding the precision of estimates of their numbers was warranted (Choat and Ayling,
1987). These species were not, however, omitted from the multivariate analyses, as such
large numbers do not pose a problem for the dissimilarity measure used here. Although
there still may be some difficulty in interpreting analyses of total numbers of individuals
among sites with different species composition, as the underwater visual sampling method
may be biased in different ways for different fish species, we nevertheless considered that
it was a measured variable worth investigating, with this caveat in mind.
2.3. Binomial deviance as a dissimilarity measure
For the fish abundance data, we developed and used a new dissimilarity measure, based
on likelihood theory. Let y1k and y2k be the count for species k in transects 1 and 2,
respectively, and let nk=( y1k+y2k). We may then consider the null hypothesis that the two
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 199
transects do not differ in their composition and relative abundance of species. That is, for
each species, we expect that half of the counts (i.e., half of nk) will fall into each transect.
The binomial deviance of our observed data from this hypothesis is defined, under
likelihood theory, as
y1k y2k 1
deviance ¼ y1k log þ y2k log À ðy1k þ y2k Þlog
nk nk 2
Thus, a useful measure of ecological dissimilarity between the two transects is the sum of
these deviances across all species. To take account of the fact that different species may be
varying on different scales, a scale-invariant measure may be obtained by considering this
quantity on a per-observation basis, which is achieved by simply dividing by nk:
X1
p
y1k
Y2k
1
d1;2 ¼ y1k log þ y2k log À ðy1k þ y2k Þlog
k¼1
nk nk nk 2
This measure, which we will refer to as the binomial deviance dissimilarity, resembles
somewhat the measure proposed by Cao et al. (1997) as an improvement on the Bray-
Curtis measure, but has the added advantage of being based in likelihood theory.
The null hypothesis from which the binomial deviance dissimilarity is derived could
also be tested using the chi-square test statistic for the equality of expected counts in the
two transects for all species. In this context, the chi-square test statistic is equivalent to the
dissimilarity measure known as the coefficient of divergence (Clark, 1952). Hence, the
binomial deviance and the coefficient of divergence should have similar properties. The
coefficient of divergence, along with the Canberra metric and the Bray-Curtis measure,
were recommended by Gower and Legendre (1986) for use with species abundance data.
For further details on the properties of these measures, see Legendre and Legendre, (1998,
e.g., p. 298).
2.4. Ordination
Metric multi-dimensional scaling (MDS) (principal coordinate analysis) on the basis of
the binomial deviance dissimilarity measure was used as an unconstrained ordination
method to visualize multivariate patterns. In addition, canonical analysis of principal
coordinates (CAP, Anderson and Willis, 2003; Anderson and Robinson, 2003) was used as
a constrained ordination procedure to visualize patterns by reference to particular
hypotheses. The CAP analyses were only done on appropriate terms found to be
significant using NPMANOVA. These analyses were also done using specialised software
written in FORTRAN.
2.5. Univariate analyses
Frequencies of occurrence were examined for each fish species across the relevant
factors of interest in the study (Appendix A). Six species of fish were abundant enough to
be analysed using univariate analysis of variance (ANOVA): Parika scaber (leather
200 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
jacket), Cheilodactylus spectabilis (red moki), Scorpis lineolatus (sweep), Upeneichthys
lineatus (goatfish), N. celidotus (spotty) and Chromis dispilus (demoiselle). The total
number of fish and the total number of species per transect were also analysed using
univariate ANOVA. Due to the predominance of zeros, extremely large variability at the
level of individual transects and the highly skewed nature of these data, normality was not
a reasonable assumption for any of these single variables except for the total number of
species. Thus, in each case, all tests were done using a permutation procedure (with 4999
permutations of appropriate units), as described for the NPMANOVA above. These tests
using permutations will be sensitive to differences in dispersion, so results should be
interpreted with caution in this regard. Where appropriate, significant terms were
investigated with a posteriori pair-wise comparisons, which also used 4999 random
permutations to obtain P-values. These univariate tests were done using the same software
as that used for the multivariate tests, but with only one variable and (appropriately) basing
the analysis on Euclidean distance. The F-ratios used for tests done in this way are
equivalent to those of traditional ANOVA (Anderson, 2001a), although the P-values are
not obtained using traditional tables.
3. Results
3.1. Effects of depth and density of kelp
The two habitats sampled did differ in their spatial distributions with respect to depth
(Fig. 2). The median depth of transects in kelp habitat was 13.5 m (inter-quartile
range = 11.5 –15 m), while for barrens habitat the median depth was 6.7 m (inter-quartile
Fig. 2. Boxplots of the depth of transects (in metres) for each habitat: barrens (n=316) and kelp (n=318).
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 201
Table 1
Relationships between the total number of fish species (tot sp) and the total number of fish (tot fish) vs. depth,
average kelp density (kelp av) or standard deviation of kelp density (kelp sd)
n r2 F P
Year 1
tot sp vs. depth 309 0.0043 1.325 0.2506
tot fish vs. depth 309 <0.0001 0.003 0.9538
tot sp vs. kelp av 159 0.0711 12.022 0.0007
tot fish vs. kelp av 159 0.0167 2.700 0.1043
tot sp vs. kelp sd 159 0.0774 13.171 0.0004
tot fish vs. kelp sd 159 0.0349 5.673 0.0184
Year 2
tot sp vs. depth 325 0.0057 1.869 0.1726
tot fish vs. depth 325 0.0157 2.691 0.0238
The total number of fish was transformed to ln(x+1) before analysis. Note that kelp densities were not recorded in
year 2.
range = 5.35– 8 m). There was no significant relationship, however, between the total
number of species and depth for either year 1 or for year 2 (Table 1). The total number of
fish (transformed to ln(x+1)) was significantly and positively related to depth for year 2
only (Table 1). There was also a significant relationship between the depth of transects and
the multivariate fish assemblages (Table 2). Depth only explained 3.5% of the variation in
the multivariate assemblage structure and only 1.6% of the variation in total numbers of
fish. Furthermore, effects of habitat (kelp vs. barrens) were significant over and above
effects of depth (Table 2). Thus, although some effects of habitat (described in greater
detail below) may be attributable to differences in depth, this analysis shows there were
significant effects of habitat on fish assemblages (e.g., structural or other differences) that
were unrelated to depth. Depth was not treated as a covariable in subsequent analyses (e.g.,
Table 3 below), as we wished to examine overall effects of habitat, including those aspects
that may have been related to depth.
Within kelp forest habitats for year 1, fish assemblages were also significantly related to
the average kelp density per transect (pseudo F1, 157=4.788, P=0.0008, 4999 permuta-
tions), but were not related to the variation (standard deviation) in kelp density along
transects (pseudo F1, 157 = 1.214, P = 0.3134, 4999 permutations). The total number of
Table 2
Sequential non-parametric multivariate multiple regression showing the relationship between multivariate fish
species abundance data (based on the binomial deviance dissimilarity measure) and depth, followed by the effect
of habitat, taking depth into account as a covariable
Source df SS MS F P
Depth 1 32.660 32.660 22.942 0.0002
Habitat/depth 1 10.050 10.050 7.128 0.0078
Residual 631 889.665 1.410
Total 633 932.375
202 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 3
Non-parametric multivariate analysis of variance of 46 fish species abundance variables, based on the binomial
deviance dissimilarity measure
Source df SS MS F P Denominator Permutable
MS units
Year (Ye) 1 8.280 8.280 0.979 0.4608 YeÂLo 8 YeÂLo units
Location = (Lo) 3 87.368 29.123 5.765 0.0002 Si(Lo) 16 Si(Lo) units
Site(Location) = 12 60.620 5.052 4.388 0.0002 Res 640 observation
Si(Lo) units
Habitat = Ha 1 25.708 25.708 2.264 0.0650 LoÂHa 8 LoÂHa units
YeÂLo 3 25.384 8.461 2.975 0.0012 YeÂSi(Lo) 32 YeÂSi(Lo)
units
YeÂSi(Lo) 12 34.130 2.844 2.470 0.0002 Res 640 observation
units
YeÂHa 1 4.117 4.117 1.598 0.1948 YeÂLoÂHa 16 YeÂLoÂHa
units
LoÂHa 3 34.065 11.355 3.684 0.0002 Si(Lo)ÂHa 32 Si(Lo)ÂHa
units
Si(Lo)ÂHa 12 36.983 3.082 2.677 0.0002 Res 640 observation
units
YeÂLoÂHa 3 7.728 2.576 1.256 0.2248 YeÂSi(Lo)ÂHa 64 YeÂSi(Lo)Â
Ha units
YeÂSi(Lo)ÂHa 12 24.609 2.051 1.781 0.0004 Res 640 observation
units
Residual 576 663.153 1.151
Total 639 1012.145
P-values were obtained using 4999 permutations of given permutable units for each term or using 4999 Monte
Carlo samples from the asymptotic permutation distribution (given in italics) when there were few possible
permutations.
species was significantly positively related to the average and to the standard deviation in
kelp density per transect, while the total number of fish (transformed to ln(x+1)) was
significantly positively related only to the standard deviation in kelp density (Table 1).
However, these relationships, although statistically significant, were very weak, with a
large amount of scatter and very low values of r2 (Table 1). This indicates that kelp density
and depth, while apparently contributing to explain small amounts of observed variation in
fish assemblages, do not provide the means to generate predictions at the scale of
individual transects.
3.2. Measured variation
Non-parametric MANOVA on the fish data showed that there was significant small-
scale variability in the fish assemblages from site to site and year to year in different
habitats (i.e., a significant YearÂSite(Location)ÂHabitat interaction, Table 3). This
significant small-scale variation was reflected in the analysis of several of the
relatively abundant individual species as well, namely: N. celidotus, S. lineolatus
and C. dispilus (see Tables 6 and 8 below). The greatest multivariate variability
occurred at the scale of individual transects (i.e., residual, Table 4). The next-greatest
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 203
Table 4
Estimated multivariate pseudo variance components for each term in the model based on sums of squared
dissimilarities (binomial deviance) for 46 fish species and the analogous univariate ANOVA estimators using
mean squares (i.e., using multivariate mean squares in Table 1), ordered from largest to smallest
Source MS Variance components
Residual 1.1513 1.1513
Lo 29.1227 0.1504
LoÂHa 11.3551 0.1034
Si(Lo) 5.0516 0.0975
Si(Lo)ÂHa 3.0819 0.0965
YeÂSi(Lo)ÂHa 2.0507 0.0899
YeÂSi(Lo) 2.8441 0.0846
Ha 25.7084 0.0449
YexLo 8.4614 0.0702
YeÂLoÂHa 2.5760 0.0131
YeÂHa 4.1171 0.0096
Ye 8.2796 0.0000
All components are random variance components except for habitat, year and YeÂHa which are sums of squared
fixed effects.
measured source of variation was that across locations, followed by the random
interaction of location with habitat, and then the variation from site to site within
locations. Considerable variation was also contributed by the random interaction of
sites with habitats (Table 4). Other measured sources of variation included interactions
of spatial effects with years; however, overall changes in assemblages from year to
year were not strong compared to the other effects. In fact, the estimated variation due
to years was slightly negative (À0.0006) using the ANOVA estimator and thus was
estimated as zero (Table 4).
3.3. Effects of habitat
The effects of habitat (kelp forest vs. barrens) on fish assemblages varied significantly
across locations (Table 3). In particular, the effects of habitat appeared to be strongest at
Hahei and Home Point. This is shown by a greater separation of the points representing
assemblages in kelp vs. barrens assemblages for these two locations, compared to Berghan
Point and Leigh in the unconstrained ordinations: the two-factor metric MDS plot (Fig. 3)
and in each of the one-factor metric MDS plots (Fig. 4, left-hand side). In addition, the
CAP analysis showed larger canonical correlation coefficients for Hahei and Home Point
than for Leigh or Berghan Point (Table 5, Fig. 4, right-hand side). Furthermore, the
allocation success showed that fish assemblages from the two habitats at either Leigh or
Berghan Point were not as distinct as that seen for the other locations. For example, Leigh
had only 69% allocation success compared to Hahei, which had 94% success (Table 5).
Pair-wise comparisons showed that there were significant differences in assemblages of
fishes between barrens and kelp habitats at either Hahei or Home Point, but not at Berghan
Point or Leigh (Table 5). It is interesting to note that the CAP analysis for Leigh shows a
fairly distinct separation due to habitat, although its allocation success was poor and there
was no significant effect detected by NPMANOVA. This is essentially due to the selection
204 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Fig. 3. Two-factor metric MDS plot of the pooled fish assemblages for each site on the basis of the binomial
deviance dissimilarity measure, showing the factors of location and habitat. There are eight points for each
combination of the factors, which correspond to the four sites in each of 2 years.
effects occurring in high-dimensional space, which can cause canonical plots to paint a
‘‘rosy’’ picture of the results. It further emphasizes the importance of using a technique
like the leave-one-out method (e.g., Lachenbruch and Mickey, 1968) for assessing the
distinctness of groups and independent tests rather than relying solely on the canonical
plot.
There were no consistent effects of habitat on the total number of fish; these effects
varied from site to site and from year to year (Table 6a, Fig. 5a). Despite the
significant three-way interaction of YearÂLocationÂHabitat (Table 6b), there were no
statistically significant pair-wise differences detected in the total number of species in
kelp forest vs. barrens in either year at any location (pair-wise comparisons, P > 0.05,
Fig. 5b). In contrast, average numbers of P. scaber (leather jackets) were significantly
higher in kelp forests than in barrens habitats for both years at the two northern
locations of Berghan Point and Home Point (Table 6c, Fig. 5c, P < 0.05). Their
frequency of occurrence was also greater in kelp forests, as was that of C. dispilus
(demoiselle), T. novaezelandiae (jack mackerel), Nemadactylus douglasii (porae),
Bodianus unimaculatus (pigfish), Odax pullus (butterfish) and Pseudolabrus miles
Fig. 4. Unconstrained metric MDS plots (left) and constrained CAP plots (right) done separately at each location
(rows) on the basis of the binomial deviance dissimilarity measure, in each case comparing fish assemblages in
two different habitats: barrens vs. kelp. There are eight points for each combination of the factors, which
correspond to the four sites in each of 2 years.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 205
206 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 5
Results of CAP analyses examining effects of habitat within each location for 46 species of fish on the basis of the
binomial deviance dissimilarity measure
Location m %Var Allocation success (%) d2 P
Barrens Kelp Total
Berghan Point 4 78.56 75 75 75 0.607 0.0564
Home Point 5 84.81 88 75 81 0.748 0.0018
Leigh 6 86.97 63 75 69 0.690 0.1080
Hahei 4 86.89 88 100 94 0.824 0.0006
m = the number of principal coordinate (PCO) axes used in the CAP procedure, %Var = the percentage of the total
variance explained by the first m PCO axes, Allocation success = the percentage of points correctly allocated into
each group, d2 = the squared canonical correlation. P-values given are the results of pair-wise comparisons of
assemblages in kelp forest vs. barrens habitat at each location, using NPMANOVA, with 4999 permutations of
individual sites as units (i.e., permuting each set of n = 10 transects together as a unit).
(scarlet wrasse) (Appendix A). Correlations of species with canonical axes also
suggested that Zeus faber (john dory) and Seriola lalandi (kingfish) were associated
with kelp forest habitats (Table 7). Although average numbers of C. dispilus also
tended to be greater in kelp forest habitat (Fig. 5d), this pattern was not statistically
significant (Table 6d).
Fish that occurred more frequently in barrens habitats were N. celidotus (spotty),
Notolabrus fucicola (banded wrasse), Girella tricuspidata (parore), Coris sandageri
Table 6
Results of univariate ANOVAs on selected variables and species
Source df (a) Total no. fish (b) Total no. species (c) Parika scaber (d) Chromis dispilus
MS F P MS F P MS F P MS F P
Ye 1 8917.0 4.779 0.1114 3.025 0.070 0.8108 2.336 0.892 0.4206 4944.4 4.565 0.1164
Lo 3 13056.2 4.065 0.0370 98.173 5.342 0.0170 15.135 2.789 0.0888 9972.5 7.064 0.0088
Si(Lo) 12 3211.7 4.916 0.0002 18.377 4.853 0.0002 5.427 2.674 0.0014 1411.7 4.575 0.0002
Ha 1 17.2 0.010 0.9332 32.400 0.388 0.5752 43.403 2.908 0.1788 3048.9 4.761 0.1084
YeÂLo 3 1866.0 1.934 0.1824 43.050 9.156 0.0026 2.618 1.329 0.3076 1083.2 1.516 0.2600
YeÂSi(Lo) 12 965.1 1.477 0.1268 4.702 1.242 0.2538 1.970 0.970 0.4858 714.4 2.315 0.0068
YeÂHa 1 1196.5 0.586 0.5122 12.656 0.762 0.4538 9.344 8.959 0.0572 3032.5 5.714 0.0948
LoÂHa 3 1794.5 0.564 0.6510 83.508 7.402 0.0032 14.924 6.991 0.0036 640.4 0.738 0.5406
Si(Lo)ÂHa 12 3184.4 4.874 0.0002 11.281 2.979 0.0010 2.135 1.052 0.4020 867.2 2.810 0.0022
YeÂLoÂHa 3 2043.5 1.335 0.3046 16.606 4.516 0.0278 1.043 0.900 0.4670 530.7 0.523 0.6790
YeÂSi(Lo)Â 12 1530.7 2.343 0.0062 3.677 0.971 0.4762 1.159 0.571 0.8830 1013.9 3.286 0.0002
Ha
Residual 576 653.4 3.787 2.030 308.6
Total 639
P-values were obtained by 4999 permutations of appropriate units, as shown in Table 3 for each term in the
analysis.
P-values in italics were obtained using 4999 Monte Carlo samples from the asymptotic permutation distribution.
Mean squares used in the denominator for each test are also shown in Table 3.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 207
Fig. 5. Means (F1 S.E.) for (a) the total number of fish, (b) the total number of species, (c) the number of P.
scaber (leather jacket) and (d) the number of C. dispilus (demoiselle) in each combination of
YearÂLocationÂHabitat (n = 4 sites per combination of levels and a site consisted of counts summed across
10 transects).
(Sandager’s wrasse), Chironemus marmoratus (hiwihiwi), Parma alboscapularis
(black angelfish), Scorpis violaceus (blue maomao) and Kyphosus sydneyanus (silver
drummer) (Appendix A). All of these fish had strong correlations with canonical
axes, as did Pempheris adspersus (big eye) (Table 7). Average numbers of N.
celidotus also tended to be greater in barrens habitats (Fig. 6a), although this trend
was only statistically significant at Hahei ( P < 0.05), but not at the other locations
(Table 8a).
Some of the common species did not show any strong consistent patterns or differences
across the two habitats, including U. lineatus (goatfish, Fig. 6b, Table 8b), S. lineolatus
(sweep, Fig. 6c, Table 8c) and C. spectabilis (red moki), which was significantly more
208 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 7
Correlations of individual species with the canonical axis for habitat for each of the four locations, as shown in
Fig. 4 (right-hand side)
Name Berghan Home Leigh Hahei Average
Positive correlation (kelp)
Bodianus unimaculatus 0.489 0.639 – 0.232 0.453
Seriola lalandi 0.499 – 0.218 0.561 0.426
Parika scaber 0.681 0.856 0.384 À0.226 0.424
Zeus faber 0.356 0.425 0.475 0.417 0.418
Nemadactylus douglasii 0.491 0.776 0.333 À0.245 0.338
Trachurus novaezelandiae 0.420 0.305 0.250 0.060 0.259
Pseudolabrus miles – 0.270 0.075 0.398 0.247
Chromis dispilus 0.435 0.286 À0.269 0.237 0.172
Odax pullus 0.104 0.456 0.248 À0.172 0.159
Negative correlation (barrens)
Parma alboscapularis À0.409 À0.495 – À0.539 À0.481
Chironemus marmoratus À0.552 À0.023 À0.681 À0.591 À0.462
Coris sandageri À0.352 À0.303 À0.438 À0.551 À0.411
Notolabrus celidotus À0.340 À0.083 À0.380 À0.694 À0.374
Girella tricuspidata À0.414 À0.499 À0.056 À0.442 À0.353
Pempheris adspersus À0.124 À0.102 À0.569 À0.261 À0.264
Notolabrus fucicola À0.119 0.062 À0.028 À0.822 À0.227
Kyphosus sydneyanus À0.125 À0.358 0.218 À0.488 À0.188
Scorpis violaceus 0.347 À0.415 À0.090 À0.456 À0.153
A positive correlation indicates species associated with kelp habitat, while a negative correlation indicates species
associated with barrens habitat. Species are given in decreasing order of the absolute value of their average for the
correlation across the four locations. Dashed lines indicate a species did not occur at that location. Species that
occurred in fewer than 7 transects (out of a total of 640) were not included.
abundant in barrens habitats but only in year 2 and only at Hahei (Fig. 6d, Table 8d, pair-
wise comparisons, P < 0.05).
3.4. Effects of locations
Although locations were originally chosen randomly and so treated as a random
factor in the analyses, it was nevertheless of interest to compare the 4 locations in
terms of the fish assemblages found there, for biogeographic reasons and to compare
results with those of previous studies. Due to its interaction with habitat (Table 3), the
potential differences among locations were considered separately for each habitat. For
each habitat, unconstrained ordinations did not separate assemblages from different
locations very clearly (non-metric MDS plots, Fig. 7, right-hand side), while con-
strained (CAP) ordinations appeared to successfully separate Leigh, Hahei and Home
Point, with assemblages from Berghan Point being less distinct (Fig. 7, right-hand
side). The separation in multivariate space of assemblages from different locations was
slightly more successful for barrens habitats than for kelp forests (Table 9). For barrens
habitats, fish assemblages from each location differed significantly from all other
locations, except for Berghan Point and Home Point, which did not differ significantly
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 209
Fig. 6. Means (F1 S.E.) for (a) the number of N. celidotus (spotty), (b) the number of U. lineatus (goatfish),
(c) the number of S. lineolatus (sweep) and (d) the number of C. spectabilis (red moki) in each combination of
YearÂLocationÂHabitat (n = 4 sites per combination of levels and a site consisted of counts summed across 10
transects).
(Table 10). For kelp forest habitats, fish assemblages did not differ between Berghan
Point and Hahei, but all other comparisons among locations were statistically
significant (Table 10).
The average total number of fish observed per transect was significantly greater at the
two northern locations, Berghan Point and Home Point, than at either Leigh or Hahei (pair-
wise comparisons, P<0.05, Table 6a, Fig. 5a). The greatest average number of species was
found at Hahei for barrens habitats and at Home Point for kelp forests (Table 6b, Fig. 5b).
In year 1, Home Point had a significantly greater average number of species than Leigh or
Hahei, for kelp or barrens habitats. In year 2, this was only true for kelp habitats: in barrens
210 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 8
Results of univariate ANOVAs on selected species
Source df (a) Notolabrus (b) Upeneichthys (c) Scorpis (d) Cheilodactylus
celidotus lineatus lineolatus spectabilis
MS F P MS F P MS F P MS F P
Ye 1 0.046 0.003 0.9584 17.150 0.411 0.5706 228.1 0.398 0.5770 0.693 0.277 0.6396
Lo 3 28.859 6.164 0.0114 97.709 5.018 0.0038 249.9 0.675 0.6062 6.076 3.401 0.0506
Si(Lo) 12 4.682 2.122 0.0126 19.471 2.397 0.0044 370.4 2.275 0.0072 1.787 2.248 0.0104
Ha 1 71.939 3.579 0.1430 15.242 2.524 0.2014 152.0 0.224 0.6660 2.132 0.358 0.5912
YeÂLo 3 15.374 5.017 0.0204 41.756 1.506 0.2226 573.7 1.733 0.2080 2.499 2.860 0.0786
YeÂSi 12 3.064 1.388 0.1678 27.726 3.414 0.0002 331.1 2.034 0.0142 0.874 1.099 0.3466
(Lo)
YeÂHa 1 0.858 0.105 0.8190 4.350 2.096 0.2594 232.9 0.684 0.4786 2.634 0.753 0.4688
LoÂHa 3 20.103 3.849 0.0334 6.040 0.963 0.4586 678.0 1.391 0.2984 5.949 5.654 0.0126
Si(Lo)Â 12 5.223 2.366 0.0064 6.270 0.772 0.7356 487.3 2.994 0.0008 1.052 1.324 0.2068
Ha
YeÂLo 3 8.201 2.050 0.1578 2.076 0.402 0.7866 340.6 0.984 0.4448 3.497 6.829 0.0048
Ha
YeÂSi 12 4.001 1.813 0.0436 5.166 0.636 0.8532 346.3 2.127 0.0068 0.512 0.644 0.8154
(Lo)Ha
Residual 576 2.207 8.123 162.8 0.795
Total 639
P-values were obtained by 4999 permutations of appropriate units, as shown in Table 3 for each term in the
analysis.
P-values in italics were obtained using 4999 Monte Carlo samples from the asymptotic permutation distribution.
Mean squares used in the denominator for each test are also shown in Table 3.
habitats, Hahei had a significantly greater average number of fish species than either
Berghan Point or Leigh (Fig. 5b, pair-wise comparisons).
Significantly greater numbers of P. scaber were found at Hahei than at other
locations for barrens habitats, while for kelp forests, there were significantly greater
average numbers at Home Point compared to the other locations (Fig. 5c, Table 6c).
There were no significant differences among locations in numbers of S. lineolatus, but
C. dispilus was significantly more abundant at the two northern locations: Berghan
Point and Home Point, compared to Leigh and Hahei (Fig. 5d, Table 6d, pair-wise
comparisons, P<0.05). N. celidotus was significantly more abundant at Hahei than at
the other locations in barrens habitats, while there were no significant differences in its
average abundances at different locations for kelp habitats (Table 8a, Fig. 6a, pair-wise
comparisons). U. lineatus was significantly more abundant at Leigh than at the other
locations (Table 8b, Fig. 6b). The only statistically significant differences among
locations for C. spectabilis occurred in year 2 and in barrens habitat only, where
Hahei had significantly greater numbers, on average, than either Leigh or Berghan
Point (Fig. 6d, Table 8d).
Some species were more frequently observed at Home Point than at any other locations,
including P. scaber, P. alboscapularis and N. douglasii (Appendix A). Others had greater
frequencies of occurrence at the two northern locations (Berghan Point and Home Point)
than at the southern locations (Leigh and Hahei), such as C. dispilus, C. sandageri and B.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 211
Fig. 7. Unconstrained metric MDS plots (left) and constrained CAP plots (right) done separately for each habitat
(rows) on the basis of the binomial deviance dissimilarity measure, in each case comparing fish assemblages
among four different locations: Berghan Point, Home Point, Leigh and Hahei. There are eight points for each
combination of the factors, which correspond to the four sites in each of 2 years.
unimaculatus (Appendix A). Alternatively, some species were observed significantly more
frequently at Home Point and Hahei than at any other location: C. spectabilis, G.
tricuspidata, C. marmoratus and Aplodactylus arctidens (marblefish) (Appendix A).
Table 9
Results of CAP analyses examining effects of location within each habitat for 46 species of fish on the basis of the
binomial deviance dissimilarity measure
Habitat m %Var Allocation success (%) d2
B Ho L Ha Total
Barrens 8 91.91 88 63 75 88 78 0.734
Kelp 8 94.56 63 88 50 63 66 0.697
Headings are as defined for Table 5, with B=Berghan Point, Ho=Home Point, L=Leigh, Ha=Hahei.
212 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Table 10
P-values for NPMANOVA pair-wise comparisons among locations for each Habitat, using 4999 permutations
Habitat Pair-wise P
comparison
Barrens B vs. Ho 0.0618
B vs. L 0.0228
B vs. Ha 0.0106
Ho vs. L 0.0002
Ho vs. Ha 0.0014
L vs. Ha 0.0008
Kelp B vs. Ho 0.0228
B vs. L 0.0050
B vs. Ha 0.0608
Ho vs. L 0.0002
Ho vs. Ha 0.0004
L vs. Ha 0.0006
In each case, the n=10 transects within a site were permuted together as a unit. For each test, there were eight such
units per group (four sites observed per LoÂHa combination in each of 2 years). No corrections have been made
for multiple tests. B=Berghan Point, Ho=Home Point, L=Leigh, Ha=Hahei.
Finally, U. lineatus occurred more frequently at Leigh than anywhere else, and N. fucicola
and T. novaezelandiae occurred more frequently at Hahei than anywhere else (Appendix A).
4. Discussion
There were significant differences in the multivariate structure of fish assemblages
between kelp forest and barrens habitats. These were not independent of species
identity. Differences between habitats included greater frequencies and/or abundances
of P. scaber, C. dispilus, T. novaezelandiae, N. douglasii, B. unimaculatus, O. pullus
and P. miles in kelp forest habitat, while N. celidotus, N. fucicola, G. tricuspidata, C.
sandageri, C. marmoratus, P. alboscapularis, S. violaceus and K. sydneyanus were
more frequent or abundant in barrens habitat. Some of the more common species,
including U. lineatus, S. lineolatus and C. spectabilis showed no strong consistent
effects of these two differing habitats on their distributions. These results differ
somewhat from the results of Choat and Ayling (1987), who focused instead on the
size classes and feeding groups, rather than on individual species, in their description.
The present work examines species’ overall patterns with a broad brush, while the
work by Choat and Ayling (1987) recognized the importance of ontogenetic changes
in diet and habitat use by individual species (e.g., Clements and Choat, 1993; Moran
and Clements, 2002).
One of the most likely reasons for differences in fish assemblages between habitats
is due to depth. Kelp forest and barrens habitats in northeastern New Zealand differ in
their depth distributions (Fig. 2 and see Schiel, 1990). Previous studies have shown
how the numbers and diversity of fish are affected by depth in temperate New
Zealand (e.g., Meekan and Choat, 1997; Brook, 2002) and in tropical coral reefs (e.g.,
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 213
Roberts and Ormond, 1987; Friedlander and Parrish, 1998). We found that depth had
a significant effect on fish assemblage composition (Table 2), although richness and
the total numbers of fish were not significantly related to depth (Table 1). It was clear
in the present study, however, that effects of habitat were not limited to effects of
depth alone (Table 2), suggesting that other aspects of the habitat were important in
structuring assemblages.
It may well be that fundamental biological differences in diet and patterns of foraging
are driving differences in relative frequencies or abundances of species in different
habitats, as suggested by Choat and Ayling (1987). For some species (e.g., O. pullus),
the kelp E. radiata provides a source of food, although Clements and Choat (1993)
suggested that the brown alga Carpophyllum spp. may be preferred by this species. In
contrast, although K. sydneyanus is a herbivore that consumes Ecklonia (Moran and
Clements, 2002), it was found more frequently in barrens habitat in the present study
(Appendix A).
There are many other possible mechanisms that might also explain differences in fish
assemblages in different habitats. There are natural differences in habitat complexity
between kelp and barrens habitats, with kelp forests providing a more complex three-
dimensional biogenic structure. As such, the kelp may provide a refuge from predation,
particularly for juvenile stages (e.g., P. scaber). Also, within kelp forests, we found that
diversity and total abundance of fish were each positively related to variation in kelp
density (Table 1), indicating that the structure and heterogeneity of the habitat may play
a role. Habitat complexity can modify the effects of predation on fish by providing
refuges (e.g., Connell and Jones, 1991; Hixon and Beets, 1993; Caley and St. John,
1996; Beukers and Jones, 1997; Tupper and Boutilier, 1997; Steele, 1999). Alternatively,
differences in kelp density can result in differences in understorey algae and invertebrate
species, which fish may respond to positively or negatively (e.g., Carr, 1989). Connell
(2002) and Wellenreuther and Connell (2002) have demonstrated such effects of habitat-
driven prey availability on spatial patterns of a temperate reef fish (Cheilodactylus
nigripes) from local through to broad scales (hundreds of kilometers). Habitat may also
affect the frequency or intensity of inter-or intra-specific behavioural interactions (e.g.,
Levin et al., 2000).
It is important not to infer too much concerning processes that may be underlying
observed patterns in the results reported here. What we have given is a ‘‘snapshot’’ of
fish we happened to see when and where we saw them. First of all, not all potential
habitats for these fish have been included in this study (e.g., sponge gardens, shallow
areas with dense algal stands of Carpophyllum spp. or other algae, etc.). Thus, for
example, observing increased frequencies of some species in Ecklonia forests should not
be interpreted to mean that this is their ‘‘preferred’’ habitat. Furthermore, we expect
individual species will occur in different habitats at different stages in their life history
or during different behavioural stages (e.g., nesting, feeding, etc.), indeed even during
different stages of the tide or time of day. Quantitative observations we have given here
will have some biological basis, but this will require further and more detailed species-
specific studies.
One important result obtained here is that effects of habitat, although relatively
consistent for the 2 years of observation, did not occur at all locations. In fact,
214 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
multivariate effects of habitat were marginal at Berghan Point ( P = 0.0564) and were not
statistically significant ( P > 0.10) for Leigh (Table 5). This could be one reason that
Choat and Ayling (1987) concluded that habitat effects were independent of species
identity, as many of their survey sites were around the Leigh area. We did find, however,
that where effects of habitat were evident, the direction (i.e., the nature) of multivariate
effects was similar (e.g., Fig. 3). This allowed us to make some generalisations about the
nature of habitat differences, in terms of compositional changes, where they did occur.
We must be cautious, however, in invoking any generalisations of processes that might
have produced observed patterns. For example, Fowler-Walker and Connell (2002) have
suggested that top-down processes are not as important in shaping shallow subtidal
benthic coastal assemblages in South Australia as in New South Wales, because South
Australia lacks the predominance of herbivorous grazers found in regions of New South
Wales. Similarly, variations in effects of habitat among locations within New Zealand
could be due to variation in the abundance or activity of invertebrate grazers and the
relative strength and interactions of bottom-up vs. top-down processes (e.g., Menge,
1992).
Differences in fish assemblages among locations provided interesting biogeographic
information that supports results from previous studies and provided some new insights.
First, there were, on average, greater numbers of species at Home Point than at any
other location, followed by Berghan Point and Hahei, with Leigh, on average, having
the most depauperate fish communities. This is consistent with previous studies of North
Island New Zealand fish fauna, where higher diversity at the Poor Knights Islands, the
Karikari peninsula and Cape Brett was hypothesized to be a consequence of the
influences of the East Auckland Current (e.g., Denham et al., 1984; Choat and Ayling,
1987; Brook, 2002), which does not appear to have much of an influence further south
(i.e., at Leigh). However, in this study we found that, for kelp forest habitats, there was
no significant difference between fish assemblages at Hahei (the southern-most location)
and those at Berghan Point (the northern-most location). So, the influences of the East
Auckland Current appear to extend, at least potentially, to areas as far south as Hahei,
where species such as B. unimaculatus (pigfish), Suezichthys aylingi (crimson cleaner-
fish) and Suezichthys arquatus (rainbowfish) were all recorded in the present study. In
contrast, Leigh may exist as a kind of oceanographic ‘‘backwater’’ of biodiversity for
fish, potentially because of increases in turbidity and/or decreases in relative exposure as
one moves southward into the Hauraki Gulf (Grace, 1983). It is also possible that there
is greater fishing pressure around Leigh, which is not as isolated from human
populations as the other locations in the study. Roberts et al. (1992) also found a
non-linear gradient of fish assemblage structure with latitude in the Red Sea that was
apparently caused by (or at least related to) a negative association of diversity with
turbidity.
The analysis of multivariate variability in fish assemblages at different spatial scales
revealed that the greatest variation occurred at the smallest spatial scale, between
individual transects (Table 4). This is not terribly surprising, given that the spatial
scale of individual transects is not large compared to the high mobility of many fish
species included in these surveys. This result concurs with many studies of inverte-
brates and algae in intertidal and subtidal environments, which have also often found
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 215
the greatest variability to occur at small spatial scales (e.g., Archambault and Bourget,
1996; Underwood and Chapman, 1996; Menconi et al., 1999; Fowler-Walker and
Connell, 2002). In addition, we found that some of the most abundant and ubiquitous
species in the surveys had significant variability from year to year at different sites
and in different habitats (i.e., a significant YearÂSite(Locations)ÂHabitat term for each
of C. dispilus, N. celidotus and S. lineolatus), meaning that predictability at small
scales could prove very difficult for these species. Also, the spatial distribution of
planktivores, such as C. dispilus and S. lineolatus, will likely be heavily dependent on
temporal changes in currents, seasons and tides. We sampled only over a 1-month
period in each of 2 years, but it is clear that variation at different temporal scales
(e.g., tidal, within-day, daily, monthly, etc.) can also be quite large and can pose
difficulties in the assessment of effects of habitat for fish (e.g., Labridae, Connell and
Kingsford, 1998).
Variation from site to site in multivariate assemblages was about the same size as
variability from location to location (Table 4), which concurs with the results of Choat
and Ayling (1987). Observed spatial patterns in fish assemblages were consistent in the
2 years. These results are consistent with previous studies. Jones (1984b) indicated that,
generally, variation in the densities of N. celidotus in different sites was greater than
year-to-year variation over 4 years. Doherty (1987) indicated that, although temporal
variation was detected at virtually every spatial scale found, recruitment was relatively
consistent and hence predictable at geographic and regional scales and Fowler et al.
(1992) found consistent spatial patterns of recruitment from year to year in butterflyfish
on the Great Barrier Reef. In contrast, Tolimieri (1995) found no consistent patterns in
recruitment of damselfish from year to year in the Caribbean and Sale et al. (1984)
found significant year-by-reef interactions in recruitment for nine species of fish on the
Great Barrier Reef. The number of years examined to date (only two) is still too small to
make any inferences about inter-annual variation. Observations of New Zealand’s fish
fauna are needed over longer periods of time (i.e., over several years) and at various
temporal scales to begin to understand potential spatio-temporal interactions or consis-
tencies of patterns.
The observations of patterns given here and by Choat and Ayling (1987) provide a
starting point for further investigations of effects of habitat on distributions of
temperate New Zealand fishes (e.g., Underwood et al., 2000). Experiments, including
manipulations of habitat, observations of feeding, behaviour, competition and preda-
tion, are needed for individual species and groups to further investigate the causal
mechanisms behind any observed relationships or patterns. For example, Jones
(1984a), Carr (1989), Levin and Hay (1996) and Wellenreuther and Connell (2002)
have experimentally altered the presence and/or density of algae to investigate effects
on temperate reef fish. In New Zealand, Jones (1984a) showed experimentally how
juvenile spotties, N. celidotus, recruited in greater abundance in kelp forests and were
negatively affected by decreases in kelp density. Although variation in recruitment
among sites was apparently reflected in adult densities, Jones (1984b) reported
nevertheless, on average, between 4.4 and 7.2 spotties (per 500 m2) in kelp forest
habitat sites, compared to 13.4 and 46.8 (per 500 m2) in shallow broken rock habitat
near Leigh (i.e., Table II therein). Thus, distributions of adult N. celidotus as reported in
216 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Jones (1984b) support the results found here, that is, of greater numbers, on average, in
barrens habitats.
Work along the coast of southern California and reviews of studies in other parts of the
world have emphasised the idea that temperate reef fish, particularly through their
predatory impacts on herbivores, such as urchins, can affect the distributions of kelp
forests (e.g., Cowen, 1983; Holbrook et al., 1990; Tegner and Dayton, 2000). In contrast,
previous studies in New Zealand (including the present study) have emphasised the idea
that temperate reef fishes can respond to changes in habitat, and have not found much
evidence to support the alternative idea that, instead, fish are affecting the structure and
demography of the habitat (e.g., Andrew and Choat, 1982; Choat and Ayling, 1987; Jones,
1988; Schiel, 1990). However, recent studies of the large-scale effects of marine reserves
have revealed that changes in densities of predatory fishes may have strong top-down
effects on community structure (Babcock et al., 1999; Shears and Babcock, 2002). In
particular, increases in predators, such as snapper (Pagrus auratus) and blue cod (Para-
apercis colias), inside of marine reserves resulted in increased rates of predation and
decreases in densities of urchins (E. chloroticus), which in turn has been correlated with
notable increases in the physical extent of kelp forest as opposed to barrens habitat within
reserves (Babcock et al., 1999; Shears and Babcock, 2002). Thus, perhaps not surprisingly,
it would appear that both processes occur: fish populations can affect habitat structure and
habitat, in turn, can affect fish. This recent work, combined with the present study and the
work of Choat and Ayling (1987), suggests that there may be dynamic feedback loops
between indirect effects of fish on habitat structure (especially through predation on
urchins) and direct and indirect effects of habitat structure (particularly the presence and
varying density of kelp forest) on fish communities. However, the relative strength of these
processes may be location-specific, and we must heed warnings of the potential dangers
inherent in generalisations, particularly over large spatial scales (e.g., Underwood and
Petraitis, 1993; Fowler-Walker and Connell, 2002). Clearly, more study is needed and of
more components of the fish fauna in New Zealand to examine potential mechanisms and
the relative importance of these processes.
Acknowledgements
This study could not have been done without the generous help of many diving
assistants and volunteers, including: K. Bloxham, R. Dixon, B. Doak, D. Egli, D. Feary, C.
Galpin, T. Langlois, G. Nesbit, D. Parsons, A. Rapson, M. Rixon, T. Ross-Watt, J.
Saunders, M. Schlegel, N. Tolimieri and T. Willis. We also thank Arthur Cozens, Russ
Babcock, the Leigh Marine Laboratory, Peter Keen and (what was then) the School of
Environmental and Marine Science for logistic support. Some of the data from this study
form part of an investigation into the large-scale patterns of variation in temperate reef
fishes, in collaboration with S. Connell and B. Gillanders (University of Adelaide,
Australian Research Council Discovery Project), who deserve special thanks for all of
their input into the initial experimental design. We also thank K. Clements, S. Connell and
T. Willis for valuable discussions and comments on the manuscript, and N. Shears for
assistance with drawing maps. [AU]
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 217
Appendix A
List of the 46 fish species recorded in the study and their frequencies of occurrence (a)
at each location (out of a possible 160 transects) and (b) in each habitat (out of a possible
320 transects).
Family Species Common (a) Location (b) Habitat
name Berghan Home Leigh Hahei Barrens Kelp
Aplodactylidae Aplodactylus Marblefish 1 12 1 8 14 8
arctidens
Arripidae Arripis trutta Kahawai 0 1 0 5 3 3
Carangidae Decapterus Koheru 2 5 0 2 4 5
koheru
Carangidae Pseudocaranx Trevally 2 6 1 0 6 3
dentex
Carangidae Seriola Kingfish 2 0 1 5 0 8
lalandi
Carangidae Trachurus Jack 15 3 26 42 24 62
novaezelandiae mackerel
Chironemidae Chironemus Hiwihiwi 13 22 7 20 47 15
marmoratus
Dasyatidae Dasyatis Short-tailed 1 5 3 0 3 6
brevicaudata stingray
Dasyatidae Dasyatis Long-tailed 0 3 0 0 1 2
thetidis stingray
Diodontidae Allomycterus Porcupinefish 1 2 0 0 2 1
jaculiferus
Girellidae Girella Parore 13 33 9 21 61 15
tricuspidata
Kyphosidae Kyphosus Silver 2 4 1 4 10 1
sydneyanus drummer
Labridae Anampses Elegant 1 0 0 0 1 0
elegans wrasse
Labridae Bodianus Red pigfish 8 12 0 1 4 17
unimaculatus
Labridae Coris picta Combfish 0 2 0 0 2 0
Labridae Coris Sandager’s 26 27 3 12 44 24
sandageri wrasse
Labridae Notolabrus Spotty 71 69 84 99 175 148
celidotus
Labridae Notolabrus Banded 18 24 4 52 69 29
fucicola wrasse
Labridae Pseudolabrus Orange 1 4 2 0 3 4
luculentus wrasse
Labridae Pseudolabrus Scarlet 0 5 4 9 4 14
miles wrasse
Labridae Suezichthys Rainbowfish 1 0 0 1 1 1
arquatus
(continued on next page)
218 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
Appendix A (continued)
Family Species Common (a) Location (b) Habitat
name
Berghan Home Leigh Hahei Barrens Kelp
Labridae Suezichthys Crimson 0 1 0 2 2 1
aylingi cleanerfish
Latridae Cheilodactylus Red moki 62 84 45 70 135 126
spectabilis
Latridae Latridopsis Blue moki 1 0 0 0 0 1
ciliaris
Latridae Latridopsis Copper 0 1 0 1 1 1
forsteri moki
Latridae Nemadactylus Porae 5 11 5 1 2 20
douglasii
Latridae Nemadactylus Tarakihi 0 5 0 0 0 5
macropterus
Monacanthidae Parika scaber Leatherjacket 47 83 28 59 82 135
Mugilidae Aldrichetta Yellow-eyed 0 3 1 0 1 3
forsteri mullet
Mullidae Upeneichthys Goatfish 40 41 70 26 88 89
lineatus
Muraenidae Gymnothorax Yellow moray 0 1 1 2 2 2
prasinus
Myliobatidae Myliobatis Eagle ray 2 6 2 1 5 6
tenuicaudatus
Odacidae Odax pullus Butterfish 4 7 1 8 7 13
Pempheridae Pempheris Big eye 14 9 10 4 21 16
adspersus
Pinguipedidae Parapercis Blue cod 3 5 1 2 8 3
colias
Pomacentridae Chromis Demoiselle 101 115 16 69 137 164
dispilus
Pomacentridae Chromis Yellow 1 0 0 0 0 1
fumea demoiselle
Pomacentridae Chromis Single-spot 1 0 0 1 1 1
hypsilepis demoiselle
Pomacentridae Parma Black 7 21 0 7 35 0
alboscapularis angelfish
Scorpaenidae Scorpaena Northern 0 0 0 2 2 0
cardinalis Scorpionfish
Scorpidae Scorpis Sweep 42 43 48 56 93 96
lineolatus
Scorpidae Scorpis Blue 3 18 2 3 20 6
violaceus maomao
Serranidae Caesioperca Butterfly 0 0 3 0 0 3
lepidoptera perch
Sparidae Pagrus Snapper 19 14 48 3 42 42
auratus
Trachichthyidae Optivus Slender 5 1 2 3 3 8
elongatus roughy
Zeidae Zeus faber John dory 2 1 2 2 0 7
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 219
References
Anderson, M.J., 2001. A new method for non-parametric multivariate analysis of variance. Aust. Ecol. 26,
32 – 46.
Anderson, M.J., 2001. Permutation tests for univariate or multivariate analysis of variance and regression. Can. J.
Fish. Aquat. Sci. 58, 626 – 639.
Anderson, M.J., Robinson, J., 2003. Generalised discriminant analysis based on distances. Aust. N.Z. J. Stat. 45,
301 – 318.
Anderson, M.J., ter Braak, C.J.F., 2003. Permutation tests for multi-factorial analysis of variance. J. Stat. Comput.
Simul. 73, 85 – 113.
Anderson, M.J., Willis, T.J., 2003. Canonical analysis of principal coordinates: a useful method of constrained
ordination for ecology. Ecology 84, 511 – 525.
Andrew, N.L., Choat, J.H., 1982. The influence of predation and conspecific adults on the abundance of juvenile
Evechinus chloroticus (Echinoidea: Echinodermata). Oecologia 54, 80 – 87.
Andrew, N.L., Mapstone, B.D., 1987. Sampling and the description of spatial pattern in marine ecology. Ocean-
ogr. Mar. Biol., Ann. Rev. 25, 39 – 90.
Archambault, P., Bourget, E., 1996. Scales of coastal heterogeneity and benthic intertidal species richness,
diversity and abundance. Mar. Ecol. Prog. Ser. 136, 111 – 121.
Ault, T.R., Johnson, C.R., 1998. Spatially and temporally predictable fish communities on coral reefs. Ecol.
Monogr. 68, 25 – 50.
Babcock, R.C., Kelly, S., Shears, N.T., Walker, J.W., Willis, T.J., 1999. Changes in community structure in
temperate marine reserves. Mar. Ecol. Prog. Ser. 189, 125 – 134.
Benedetti-Cecchi, L., 2001. Variability in abundance of algae and invertebrates at different spatial scales on rocky
sea shores. Mar. Ecol. Prog. Ser. 215, 79 – 92.
Beukers, J.S., Jones, G.P., 1997. Habitat complexity modifies the impact of piscivores on a coral reef fish
population. Oecologia 114, 50 – 59.
Brock, R.E., 1982. A critique of the visual census method for assessing coral reef fish populations. Bull. Mar. Sci.
32, 269 – 276.
Brook, F.J., 2002. Biogeography of near-shore reef fishes in northern New Zealand. J. R. Soc. N.Z. 32, 243 – 274.
Caley, M.J., St. John, J., 1996. Refuge availability structures assemblages of tropical reef fishes. J. Anim. Ecol.
65, 414 – 428.
Cao, Y., Williams, W.P., Bark, A.W., 1997. Similarity measure bias in river benthic Aufwuchs community
analysis. Water Environ. Res. 69 (1), 95 – 106.
Carr, M.H., 1989. Effects of macroalgal assemblages on the recruitment of temperate zone reef fishes. J. Exp.
Mar. Biol. Ecol. 126, 59 – 76.
Chesson, P., 1998. Spatial scale in the study of reef fishes: a theoretical perspective. Aust. J. Ecol. 23,
209 – 215.
Choat, J.H., Ayling, A.M., 1987. The relationship between habitat structure and fish faunas on New Zealand
reefs. J. Exp. Mar. Biol. Ecol. 110, 257 – 284.
Clark, P.J., 1952. An extension of the coefficient of divergence for use with multiple characters. Copeia 1952,
61 – 64.
Clarke, K.R., 1993. Non-parametric multivariate analyses of changes in community structure. Aust. J. Ecol. 18,
117 – 143.
Clements, K.D., Choat, J.H., 1993. Influence of season, ontogeny and tide on the diet of the temperate marine
herbivorous fish Odax pullus (Odacidae). Mar. Biol. 117, 213 – 220.
Connell, S.D., 2002. Effects of predator and prey on a foraging reef fish: implications for understanding density-
dependent growth. J. Fish Biol. 60, 1551 – 1561.
Connell, S.D., Jones, G.P., 1991. The influence of habitat complexity on postrecruitment processes in a temperate
reef fish population. J. Exp. Mar. Biol. Ecol. 151, 271 – 294.
Connell, S.D., Kingsford, M.J., 1998. Spatial, temporal and habitat-related variation in the abundance of large
predatory fish at One Tree Reef, Australia. Coral Reefs 17, 49 – 57.
Cowen, R.K., 1983. The effect of sheephead (Semicossyphus pulcher) predation on red sea urchin (Strong-
gylocentrotus franciscanus) populations: an experimental analysis. Oecologia 58, 249 – 255.
220 M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221
DeMartini, E.E., Roberts, D.A., 1990. The effects of giant kelp (Macrocystis) on the density and abundance of
fishes in a cobble-bottom kelp forest. Bull. Mar. Sci. 46, 287 – 300.
Denham, R.N., Bannister, K.M., Guthrie, K.M., Crook, F.G., 1984. Surveys of the east Auckland and east cape
currents, New Zealand. Aust. J. Mar. Freshw. Res. 35, 491 – 504.
Doherty, P.J., 1987. The replenishment of populations of coral reef fishes, recruitment surveys, and the problems
of variability manifest on multiple scales. Bull. Mar. Sci. 41, 411 – 422.
Fowler, A.J., Doherty, P.J., Williams, D.MCB., 1992. Multi-scale analysis of recruitment of a coral reef fish on
the Great Barrier Reef. Mar. Ecol. Prog. Ser. 82, 131 – 141.
Fowler-Walker, M.J., Connell, S.D., 2002. Opposing states of subtidal habitat across temperate Australia: con-
sistency and predictability in kelp canopy-benthic associations. Mar. Ecol. Prog. Ser. 240, 49 – 56.
Friedlander, A.M., Parrish, J.D., 1998. Habitat characteristics affecting fish assemblages on a Hawaiian coral reef.
J. Exp. Mar. Biol. Ecol. 224, 1 – 30.
´ ´
Garcıa-Charton, J.A., Perez-Ruzafa, A., 2001. Spatial pattern and the habitat structure of a Mediterranean reef
fish local assemblage. Mar. Biol. 138, 917 – 934.
Gower, J.C., Legendre, P., 1986. Metric and Euclidean properties of dissimilarity coefficients. J. Classif. 3, 5 – 48.
Grace, R.V., 1983. Zonation of sublittoral rocky bottom marine life and its changes from the outer to the inner
Haruaki Gulf, north-eastern New Zealand. Tane 29, 97 – 108.
Hixon, M.A., Beets, J.P., 1993. Predation, prey refuges and the structure of coral-reef fish assemblages. Ecol.
Monogr. 63, 77 – 101.
Holbrook, S.J., Schmitt, R.J., Ambrose, R.F., 1990. Biogenic habitat structure and characteristics of temperate
reef fish assemblages. Aust. J. Ecol. 15, 489 – 503.
Holbrook, S.J., Forrester, G.E., Schmitt, R.J., 2000. Spatial patterns in abundance of a damselfish reflect
availability of suitable habitat. Oecologia 122, 109 – 120.
Hurlbert, S.H., 1984. Pseudoreplication and the design of ecological field experiments. Ecol. Monogr. 54,
187 – 211.
Jones, G.P., 1984. Population ecology of the temperate reef fish Pseudolabrus celidotus Bloch and Schneider
(Pisces: Labridae): I. Factors influencing recruitment. J. Exp. Mar. Biol. Ecol. 75, 257 – 276.
Jones, G.P., 1984. Population ecology of the temperate reef fish Pseudolabrus celidotus Bloch and Schneider
(Pisces: Labridae): II. Factors influencing adult density. J. Exp. Mar. Biol. Ecol. 75, 277 – 303.
Jones, G.P., 1988. Ecology of rocky reef fish of north-eastern New Zealand: a review. N.Z. J. Mar. Freshw. Res.
22, 445 – 462.
Kotliar, N.B., Wiens, J.A., 1990. Multiple scales of patchiness and patch structure: a hierarchical framework for
the study of heterogeneity. Oikos 59, 253 – 260.
Lachenbruch, P.A., Mickey, M.R., 1968. Estimation of error rates in discriminant analysis. Technometrics 10, 1 – 11.
Legendre, P., Fortin, M.J., 1989. Spatial pattern and ecological analysis. Vegetatio 80, 107 – 138.
Legendre, P., Legendre, L., 1998. Numerical Ecology, 2nd English edition Elsevier, Amsterdam.
Levin, P.S., Hay, M.E., 1996. Responses of temperate reef fishes to alterations in algal structure and species
composition. Mar. Ecol. Prog. Ser. 134, 37 – 47.
Levin, P.S., Toimieri, N., Nicklin, M., Sale, P.F., 2000. Integrating individual behavior and population ecology:
the potential for habitat-dependent population regulation in a reef fish. Behav. Ecol. 11, 565 – 571.
Lincoln-Smith, M.P., 1989. Improving multispecies rocky reef fish censuses by counting different groups of
species using different procedures. Environ. Biol. Fishes 26, 29 – 37.
Manly, B.F.J., 1997. Randomization, Bootstrap and Monte Carlo Methods in Biology Chapman & Hall, London.
McArdle, B.H., Anderson, M.J., 2001. Fitting multivariate models to community data: a comment on distance-
based redundancy analysis. Ecology 82, 290 – 297.
McClanahan, T.R., Arthur, R., 2001. The effect of marine reserves and habitat on populations of east African
coral reef fishes. Ecol. Appl. 11, 559 – 569.
Meekan, M.G., Choat, J.H., 1997. Latitudinal variation in abundance of herbivorous fishes: a comparison of
temperate and tropical fishes. Mar. Biol. 128, 373 – 383.
Menconi, M., Benedetti-Cecchi, L., Cinelli, F., 1999. Spatial and temporal variability in the distribution of algae
and invertebrates on rocky shores in the northwest Mediterranean. J. Exp. Mar. Biol. Ecol. 233, 1 – 23.
Menge, B.A., 1992. Community regulation: under what conditions are bottom-up factors important on rocky
shores?. Ecology 73, 755 – 765.
M.J. Anderson, R.B. Millar / J. Exp. Mar. Biol. Ecol. 305 (2004) 191–221 221
Moran, D., Clements, K.D., 2002. Diet and endogenous carbohydrases in the temperate marine herbivorous fish
Kyphosus sydneyanus. J. Fish Biol. 60, 1190 – 1203.
Paulin, C., Stewart, A., Roberts, C., McMillan, P., 1989. New Zealand fish: a complete guide. National Museum
of New Zealand Miscellaneous Series, vol. 19. GP Books, Wellington, New Zealand, pp. 1 – 279.
Roberts, C.M., Ormond, R.F.G., 1987. Habitat complexity and coral reef fish diversity and abundance on Red Sea
fringing reefs. Mar. Ecol. Prog. Ser. 41, 1 – 8.
Roberts, C.M., Dawson Shepard, A.R., Ormond, R.F.G., 1992. Large-scale variation in assemblage structure of
Red Sea butterflyfishes and angelfishes. J. Biogeogr. 19, 239 – 250.
Sale, P.F., 1998. Appropriate spatial scales for studies of reef-fish ecology. Aust. J. Ecol. 23, 202 – 208.
Sale, P.F., Sharp, B.J., 1983. Corrections for bias in visual transect censuses of coral reef fishes. Coral Reefs 2,
37 – 42.
Sale, P.F., Doherty, P.J., Eckert, G.J., Douglas, W.A., Ferrell, D.J., 1984. Large scale spatial and temporal
variation in recruitment of fish populations to fish populations on coral reefs. Oecologia 64, 191 – 198.
Schiel, D.R., 1990. Macroalgal assemblages in New Zealand: structure, interactions and demography. Hydro-
biologia 192, 59 – 76.
Searle, S.R., Casella, G., McCulloch, C.E., 1992. Variance Components Wiley, Toronto.
Shears, N.T., Babcock, R.C., 2002. Marine reserves demonstrate top-down control of community structure on
temperate reefs. Oecologia 132, 131 – 142.
Snedecor, G.W., 1946. Statistical Methods, 4th edition Iowa State College Press, Ames, Iowa.
Steele, M.A., 1999. Effects of shelter and predators on reef fishes. J. Exp. Mar. Biol. Ecol. 233, 65 – 79.
Stephens Jr., J.S., Morris, P.A., Zerba, K., Love, M., 1984. Factors affecting fish diversity on a temperate reef: the
fish assemblage of Palos Verdes Point, 1974 – 1981. Environ. Biol. Fishes 11, 259 – 275.
Tegner, M.J., Dayton, P.K., 2000. Ecosystem effects of fishing in kelp forest communities. ICES J. Mar. Sci. 57,
579 – 589.
Thompson, A.A., Mapstone, B.D., 1997. Observer effects and training in underwater visual surveys of reef fishes.
Mar. Ecol. Prog. Ser. 154, 53 – 63.
Tolimieri, N., 1995. Effects of microhabitat characteristics on the settlement and recruitment of a coral reef fish at
two spatial scales. Oecologia 102, 52 – 63.
Tolimieri, N., 1998. Constrasting effects of microhabitat use on large-scale adult abundance in two families of
Caribbean reef fishes. Mar. Ecol. Prog. Ser. 167, 227 – 239.
Tolimieri, N., 1998. The relationship among microhabitat characteristics, recruitment and adult abundance in the
stoplight parrotfish, Sparisoma viride, at three spatial scales. Bull. Mar. Sci. 62, 253 – 268.
Tupper, M., Boutilier, R.G., 1997. Effects of habitat on settlement, growth, predation risk and survival of a
temperate reef fish. Mar. Ecol. Prog. Ser. 151, 225 – 236.
Underwood, A.J., Chapman, M.G., 1996. Scales of spatial patterns of distribution of intertidal invertebrates.
Oecologia 107, 212 – 224.
Underwood, A.J., Chapman, M.G., 1998. A method for analyzing spatial scales of variation in composition of
assemblages. Oecologia 117, 570 – 578.
Underwood, A.J., Petraitis, P.S., 1993. Structure of intertidal assemblages in different locations: how can local
processes be compared?. In: Ricklefs, R., Schluter, D. (Eds.), Species Diversity in Ecological Communities.
University of Chicago Press, Chicago, pp. 38 – 51.
Underwood, A.J., Chapman, M.G., Connell, S.D., 2000. Observations in ecology: you can’t make progress on
processes without understanding the patterns. J. Exp. Mar. Biol. Ecol. 250, 97 – 115.
Wellenreuther, M., Connell, S.D., 2002. Response of predators to prey abundance: separating the effects of prey
density and patch size. J. Exp. Mar. Biol. Ecol. 273, 61 – 71.
Wiens, J.A., 1989. Spatial scaling in ecology. Funct. Ecol. 3, 385 – 397.
Willis, T.J., 2001. Visual census methods underestimate density and diversity of cryptic reef fishes. J. Fish Biol.
59, 1408 – 1411.
Willis, T.J., Millar, R.B., Babcock, R.C., 2000. Detection of spatial variability in relative density of fishes:
comparison of visual census, angling, and baited underwater video. Mar. Ecol. Prog. Ser. 198, 249 – 260.
Yates, F., Cochran, W.G., 1938. The analysis of groups of experiments. J. Agric. Sci. 28, 556 – 580.