Robert M. Ewers, Stephen Thorpe, and Raphael K. Didham. 2007. Synergistic interactions between edge and area effects in a heavily fragmented landscape. Ecology 88:96–106.

Appendix A. A map of the study area, information on the sampling design, trapping method, the list of environmental and landscape variables and how they were quantified, statistical methodology, formulae for the fitted lines from the figures in the manuscript, and output tables and figures from the ordination procedures.

Supplementary Methods

Trapping method

We chose a simple, inexpensive flight interception trap (FIT) modified from the composite intercept trap of Basset (1988). The trap consisted of a rain cover and two 60 cm tall, perpendicular perspex intercept sheets above a 25 cm diameter collecting funnel and preserving jar (Fig. A2). The FITs were dug into the ground to act as a combined low-level intercept and pitfall trap. This trap type was chosen because it is suitable for capturing both ground-dwelling and canopy arthropods. We simultaneously sampled invertebrates in the forest canopy and purposefully used the same trap type at ground and canopy levels to allow direct comparisons of capture rates between the two levels (canopy data are not presented here).

Sampling design

By definition, it is not possible to sample at large distances into small fragments. Consequently, there was an unavoidable, positive correlation between fragment area and length of the edge gradient in our sampling design. This resulted in unequal sample sizes for the edge and area variables, with large fragments and small edge distances (located close to edges) having larger sample sizes than small fragments and large edge distances. To compensate for this problem, we used two statistical approaches to standardize our data for sample size (see Statistical analyses section). First, abundance estimates were converted to number of beetles captured per unit trap surface area per day. Thus, the effect of having more traps in large fragments is removed because fragment-level abundance is calculated as the mean abundance over all traps within that fragment. Second, estimates of species richness were converted from absolute richness to rarefied richness, to provide a measure that is independent of sample size.

Environmental variables

We measured 91 patch, site and landscape variables associated with the 233 sampling sites (Table A1).

Patch attributes

Fragment area (LOGAREA) was obtained from a geographic information system (GIS) analysis of the 1989 New Zealand TopoMap data set. Patch shape index (LOGSHAPE) for each fragment was calculated with the formula

SI = P/200(π.A)^0.5

where P is fragment perimeter (m) and A is fragment area (ha) (Patton 1975). Distance to edge (DISTCODE) was log₂ transformed and coded as negative or positive depending on whether the site was in forest or matrix respectively. We also included the absolute value of distance to edge (ABSDIST) based on the a priori reasoning that any species restricted to a narrow zone around forest edges (i.e., edge specialist species) would only be detected by their response to absolute distance from edge, rather than by their response to a linear gradient from forest to matrix. The interactions of fragment area with distance (AREAEDGE) and absolute distance (AREAABSD) to edge were included by multiplying their respective values together. Traps in the forest and matrix control sites were included as dummy variables (binary coding variables reflecting the presence or absence of each habitat category; FORCONT and MATCONT respectively).

Site attributes

Within a 10 × 10 m quadrat centered on each trap site we identified all canopy tree species (>10 cm dbh) and measured tree basal area (BIOMASS), litter biomass (LEAFLIT), the volume of standing dead wood (DEADST) and coarse woody debris (WOOD) with diameter > 10 cm, and visually estimated tree canopy height (CANOPY) and the percentage ground cover that fell into each of eight categories: bare ground, grass, shrub, fern, fine woody debris (branches with diameter 2 – 10 cm), litter (leaves and branches with diameter ≤ 2 cm), moss and water.

A total of 11 tree species were represented in the quadrats and tree community composition was reduced to two axes using Principal Coordinate Analysis (PCoA) on the basal area of trees (Anderson 2003). The first axis (VEGPCO1) explained 73.3 % of the variation in tree species community composition and reflected a gradient from matrix (i.e., the absence of trees) to forest species (mostly Nothofagus spp.). The second axis (VEGPCO2) explained a further 10.1 % of the variation and was correlated with a change in forest species from red beech (N. fusca) to mountain beech (N. solandri) dominated communities.

Similarly, the eight ground cover scores were reduced to three PCoA axes that explained 49.7, 22.9, and 15.4 % of the variation respectively. The first axis (QUADPCO1) described a transition from forest to pasture ground cover types, the second (QUADPCO2) a gradient in moss cover and the third (QUADPCO3) was correlated with shrub cover.

Litter samples (0.1 m²) were collected from the centre and four corners of each 10 × 10 m quadrat. Samples were dried at 80 ºC until constant weight was attained, weighed on an electronic balance to the nearest 0.01 g and the weights averaged to give one measure of litter biomass within each quadrat (LEAFLIT).

Habitat heterogeneity within the 10 × 10 m quadrat was estimated in several ways: by counting tree species richness (TREEDIV), by calculating the coefficient of variation of tree basal diameter (CVTREE), litter biomass (CVLIT), volume of dead standing trees (CVDEADST) and volume of coarse woody debris (CVWOOD), and by measuring ground cover diversity (QDIV) using the Shannon-Weiner diversity index

H’ = -Σ p_i × ln p_i

where p_i is the proportion of the quadrat covered by the i th ground cover type calculated from the eight ground cover scores.

Microclimate

Six measures of site microclimate were recorded at the 233 sites between 3 December 2002 and 24 February 2003 using Hobo microloggers. The average differential and coefficient of variation in air temperature (DELTAT, TVAR), humidity (DELTARH, RHVAR) and percent incident light (DELTAPFD, PFDVAR) were calculated by comparisons between measurements taken in the forest fragments and control values taken at two microloggers located in each of the deep forest and deep matrix controls.

Two microloggers were permanently located in each of the deep forest and deep matrix controls, and the remaining 21 available microloggers were moved between forest fragments to give a continuous 7-d period of measurement at each fragment (in a randomly allocated order). All distances within a given fragment were sampled at once in order to establish gradients in microclimate variation across the edge. Air temperature and relative humidity were recorded at 5-min intervals by internal Hobo sensors and light intensity was measured using an external light sensor designed at the University of Canterbury. Light sensors recorded a voltage signal at the micrologger which was proportional to ambient solar radiation, and each of the 25 sensors was individually calibrated to photon flux density (PFD).

For the purposes of analyses, air temperature (DELTAT) and relative humidity (DELTARH) differentials were calculated from the site value minus the matrix control (open area) value and averaged for all 5-min readings during the four hours after the solar transit (i.e., 12:45 – 16:45 hrs). This period represented the hottest part of the day when standing microclimate gradients were maximal. Averaging during this period avoided the overestimation of maximum microclimate differentials that can occur due to spurious variation in factors such as cloud cover. Variability in relative humidity (RHVAR) was measured as the coefficient of variation (CV) in 5-min readings during the same period. Variability in air temperature (TVAR) was expressed as the daily temperature range (maximum temperature within the solar transit + 4 hrs minus the minimum temperature within dawn – 4 hrs; i.e., 02:10 – 06:10 hrs).

For each site, percent incident light (DELTAPFD) was calculated as the site value divided by the matrix control value, and averaged for all 5-min readings within the four hour period ± 2 hrs of the solar transit (i.e., 10:45 – 14:45 h). In order to minimize the effect of sunfleck events and cloud cover on percent incident light, each 5-min light reading was averaged with the three earlier and three later 5-min values to give a 30-min sliding-window smoothing function. Variability in site light environment (PFDVAR) was measured as the CV of the 5-min readings during the same period.

For all six environmental variables, up to six daily mean values were obtained at each site and the average of these daily values used to compare differences in microclimate between sites (i.e., days were the replicates, not 5-min readings). Not all sites had the same number of days of microclimate recording, because weather conditions were variable and very cold conditions reduced the ability to detect microclimate gradients across edges. We conducted a sensitivity analysis of microclimate gradients to the prevailing macroclimate and excluded days in which the maximum temperature at the matrix control site did not exceed 17.5 °C (ca. 15% of days). Again, days in which more than 20% of the 5-min readings were missing within any 4-hr period were also excluded from analyses. Although 82% of sites had three or more valid days measurements, two sites did have zero valid days due to weather and animal disturbance, so values for these two sites were estimated by averaging the daily means of each variable at the two immediately adjacent sampling distances.

Landscape attributes

A GIS analysis was performed to calculate landscape variables that may influence species composition. For each fragment, the smallest distance between the fragment and the continuous forest was determined (NINADIST). We quantified landscape context for each of the 233 trap sites inside a nested set of circles with radii 256, 512, 1024, 2048, 4096, and 8192 m. For each site by radius combination we determined the number of forest fragments that were completely contained within the circle or intersected the circle perimeter (256NUM to 8192NUM), percentage of the landscape in forest (256FOR to 8192FOR), Shannon-Weiner landscape diversity (DIV256 to DIV8192), and forest connectivity (CON256 to CON8192) calculated by the formula

C = Σe^-dij × A_j

where A_j is the area of the jth neighboring forest fragment and d_ij is the distance from the trap site to the jth fragment (Steffan-Dewenter 2003).

Variables reflecting soil composition (calcium content, CALCIUM, phosphorus content, ACIDP, soil particle size, PSIZE, soil age, AGE, and soil parent material, INDUR) and macroclimate (annual and winter solar radiation, MAS and JUNES, mean annual and winter minimum temperature, MAT and TMIN, vapor pressure deficit, VPD, annual water deficit, DEFICIT, and annual water balance ratio, R2PET) were obtained from the Land Environments of New Zealand (LENZ) data set (Leathwick et al. 2003). In addition, LENZ characterises New Zealand climatic conditions and landforms at a 25 × 25 m scale in a four-level hierarchical classification (Leathwick et al. 2003), each of which was included in the analysis as dummy variables (P, K, P1, K1, E1, K11, P12, E14, K11D, E14C, E42A, and K11B).

Statistical analysis

Spatial autocorrelation

Fitted trend surfaces were used to remove spatial autocorrelation from the data and resulting analyses (Borcard et al. 1992). This approach is not as sophisticated as other spatial analytical techniques that are available (e.g., Borcard et al. 2004), and may remove some gradients in environmental variables that could be the source of variation in species abundances (Wagner and Fortin 2005), but the increased efficiency of these more sophisticated techniques is unlikely to have been realized for the widely spaced sampling locations in the present study. However, we recognize that our approach is likely to ascribe some legitimate species-environment gradients in beetle community composition to spatial autocorrelation, rather than to the underlying environmental variables which are the true driver of the observed patterns. We consider that using trend surfaces is a conservative approach to dealing with autocorrelated data, because in cases where environmental or spatial causes of a species distribution can not be differentiated, the pattern is ascribed to spatial autocorrelation and removed from the data set. Thus, only patterns in species abundances that can confidently be claimed to be distinct from spatial artefacts of the sampling design remain in the analysis.

Multicollinearity among environmental variables

We assessed multicollinearity among the 79 environmental variables using a correlation matrix (Appendix B). There were 15 pairs of variables that were highly intercorrelated (r > 0.9), and one of each pair was removed from the analysis leaving a total 64 environmental variables in the analysis.

Ordination forward selection procedure

Ordination analysis of the effects of the 64 environmental variables on beetle community composition at the 233 sites was conducted using CANOCO v4.02 (ter Braak and Smilauer 1997). To determine the set of variables that were significant drivers of community composition, we used a forward selection procedure to rank the environmental variables in order of importance (ter Braak and Verdonschot 1995). This process is analagous to forward stepwise regression, and significance was tested at each step using a Monte Carlo permutation test with 999 random permutations. The overall significance of the ordination was also tested with a Monte Carlo permutation test, using 5000 permutations.

Partitioning the drivers of edge and area effects

Both area and edge effects are composite factors that are driven by a myriad of intercorrelated environmental variables. However, many studies that report area and edge effects make little or no attempt to identify the mechanisms driving trends in community composition (e.g., Baz and Garcia-Boyero 1995, Abensperg-Traun and Smith 1999, Kotze and Samways 1999, 2001, Davies et al. 2001, Barbosa and Marquet 2002, Steffan-Dewenter 2002, Dangerfield et al. 2003). When mechanisms are investigated, the usual approach is to either break down community trends into individual species responses to area- and edge- related environmental gradients with multiple regression analyses (e.g., Didham et al. 1998, Magura 2002), or to report the correlations of environmental and habitat gradients with area and edge gradients (e.g., Pajunen et al. 1995, Abensperg-Traun et al. 1996, Didham et al. 1998, Martin and Major 2001, Brown and Freitas 2002) . The former approach provides valuable information on why individual species respond in a particular way, but falls short of describing environmental gradients that drive whole-community changes. The latter approach typically employs some form of canonical ordination analysis and gives an indication of how closely environmental gradients are aligned with community changes, but has two shortcomings: (1) when the number of environmental variables is large, the canonical coefficients become unstable and can not be used to interpret the results (ter Braak 1995), as was the case in our analysis, and (2) the coefficients and correlations are calculated between variables and the canonical axes, not the variables and the specific gradient in question. Consequently, community changes are correlated with a linear combination of all measured environmental variables (ter Braak 1986, 1995), rather than with the variable that is the subject of a specific hypothesis. Our approach to this problem was to partition the variance in the ordination that was explained by composite effects. This method involved two additional steps that follow a canonical ordination analysis. First, we extracted the portion of variation in community composition (expressed by the ordination axis scores) that could be directly explained by a specific hypothesis (in our case, area and edge gradients). Second, we partitioned that variation amongst the intercorrelated environmental gradients that are likely to drive the composite effect, using backwards stepwise regression. The final regression table shows clearly the relative contribution of the separate component drivers of the composite effect. The distinct advantage of this approach is that the interpretation of the results is more rigorous than the qualitative interpretation that normally follows ordination analyses.

Supplementary Results

Beetle abundance and species richness

Equations for the fitted lines in Fig. 1 are as follows, with A = habitat area and E = distance to edge: (a) Abundance = 8.61 + 6.30 × log₁₀(log₁₀(A)); (b) Abundance = 28.56 + (7.47 - 28.56) / (1 + 10^{((-5.04 - E) × 0.18)}); (c) forest traps only: Richness = 241.41 × (1 – e^{(-A^0.24)}); forest and matrix traps combined: Richness = 293.52 × (1 – e^{(-A^0.23)}); (d) Richness = (5.84 + 0.08 × E)(E < 0) + (5.64 – 0.15 × E)(E > 0); (e) forest traps only: Rarefied Richness = 55.34 + 0.24 × A; forest and matrix traps combined: Rarefied Richness = 67.81 + 1.48 × A; and (f) Rarefied Richness = (130.87 + 5.56 × E)(E < 0) + (115.31 – 5.82 × E)(E > 0).

Synergistic interactions between area and edge effects

Of the total pool of 893 species, 375 were categorized according to their preferred habitat by comparing abundance in the deep forest and deep matrix control sites with one-way ANOVA (Appendix C). The remaining 518 species could not be tested because they were not captured in either the deep forest or deep matrix control sites.

This measure of habitat specialization does not imply that species classified as specialists were completely restricted to their preferred habitat, but instead that their occurrence was invariably very rare in the non-preferred habitat. Species classified as forest specialists had an average abundance of 0.70 specimens captured per trap day (SE = 0.08) in the forest control sites, but just 0.02 (SE = 0.01) per trap day in the matrix controls. Similarly, matrix specialists had an average abundance of 0.80 (SE = 0.10) per trap day in the matrix controls, but just 0.01 (SE < 0.01) per trap day in the forest controls.

Equations for the fitted lines in Fig. 2 are as follows, with S = strength of edge effect and A = habitat area: (a) S = 3.91 + 2.59 × log₁₀(log₁₀(A) + 2) + 0.25 × (log₁₀(log₁₀(A) + 2)²; and (b) S = 0.067 + 0.048 × log₁₀(log₁₀(A) + 2) + 0.005 × (log₁₀(log₁₀(A) + 2)².

Equations for the fitted lines in Fig. 3 are as follows, with S = strength of area effect E = distance to edge: (a) S = 0.0216 + 0.0182 × -(log₂(│E│)) + 0.002 × -(log₂(│E│))² ; and (b) S = 0.0736 + 0.1274 × -(log₂(│E│)) + 0.0117 × -(log₂(│E│))² .

Equations for the fitted lines in Fig. 4 are as follows, with R = species richness and A = fragment area: (a) forest species: R = 1.3067 + 0.0588 × A; matrix species: R = 0.4718 + 0.0144 × A; (b) forest species: R = 1.3753 + 0.0112 × A; matrix species: R = 0.8464 – 0.0885 × A; (c) forest species: R = 1.0996 + 0.0794 ×A; matrix species: R = 0.8987 – 0.0286 ×A.

LITERATURE CITED

Abensperg-Traun, M., and G. T. Smith. 1999. How small is too small for small animals? Four terrestrial arthropod species in different-sized remnant woodlands in agricultural Western Australia. Biodiversity and Conservation 8:709–726.

Abensperg-Traun, M., G. T. Smith, G. W. Arnold, and D. E. Steven. 1996. The effects of habitat fragmentation and livestock-grazing on animal communities in remnants of gimlet Eucalyptus salubris woodland in the Western Australian wheatbelt. I. Arthropods. Journal of Applied Ecology 33:1281–1301.

Anderson, M. J. 2003. DISTLM v.2: a FORTRAN computer program to calculate distance-based multivariate analysis for a linear model. Department of Statistics, University of Auckland, Auckland, New Zealand.

Barbosa, O., and P. A. Marquet. 2002. Effects of forest fragmentation on the beetle assemblage at the relict forest of Fray Jorge, Chile. Oecologia 132:296–306.

Basset, Y. 1988. A composite interception trap for sampling arthropods in tree canopies. Journal of the Australian Entomological Society 27:213–219.

Baz, A., and A. Garcia-Boyero. 1995. The effects of forest fragmentation on butterfly communities in central Spain. Journal of Biogeography 22:129–140.

Borcard, D., P. Legendre, C. Avois-Jacquet, and H. Tuomisto. 2004. Dissecting the spatial structure of ecological data at multiple scales. Ecology 86:1826–1832.

Borcard, D., P. Legendre, and P. Drapeau. 1992. Partialling out the spatial component of ecological variation. Ecology 73:1045–1055.

Brown, K. S. J., and A. V. L. Freitas. 2002. Butterfly communities of urban forest fragments in Campinas, Sao Paulo, Brazil: structure, instability, environmental correlates, and conservation. Journal of Insect Conservation 6:217–231.

Dangerfield, J. M., A. J. Pik, D. Britton, A. Holmes, M. Gillings, I. Oliver, D. Briscoe, and A. J. Beattie. 2003. Patterns of invertebrate biodiversity across a natural edge. Austral Ecology 28:227–236.

Davies, K. F., B. A. Melbourne, and C. R. Margules. 2001. Effects of within- and between-patch processes on community dynamics in a fragmentation experiment. Ecology 82:1830–1846.

Didham, R. K., P. M. Hammond, J. H. Lawton, P. Eggleton, and N. E. Stork. 1998b. Beetle species responses to tropical forest fragmentation. Ecological Monographs 68:295–323.

Kotze, D. J., and M. J. Samways. 1999. Invertebrate conservation at the interface between the grassland matrix and natural Afromontane forest fragments. Biodiversity and Conservation 8:1339–1363.

Kotze, D. J., and M. J. Samways. 2001. No general edge effects for invertebrates at Afromontane forest/grassland ecotones. Biodiversity and Conservation 10:443–466.

Leathwick, J., G. Wilson, D. Rutledge, P. Wardle, F. Morgan, K. Johnston, M. McLeod, and R. Kirkpatrick. 2003. Land Environments of New Zealand. David Bateman Ltd., Auckland, New Zealand.

Magura, T. 2002. Carabids and forest edge: spatial pattern and edge effect. Forest Ecology and Management 157:23–37.

Martin, T. J., and R. E. Major. 2001. Changes in wolf spider (Araneae) assemblages across woodland-pasture boundaries in the central wheat-belt of New South Wales, Australia. Austral Ecology 26:264–274.

Patton, D. R. 1975. A diversity index for quantifying habitat "edge". Wildlife Society Bulletin 3:171–173.

McArdle, B. H., and M. J. Anderson. 2001. Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82:290–297.

Pajunen, T., Y. Haila, E. Halme, J. Niemelä, and P. Punttila. 1995. Ground-dwelling spiders (Arachnidae, Araneae) in fragmented old forests and surrounding managed forests in southern Finland. Ecography 18:62–72.

Steffan-Dewenter, I. 2002. Landscape context affects trap-nesting bees, wasps, and their natural enemies. Ecological Entomology 27:631–637.

Steffan-Dewenter, I. 2003. Importance of habitat area and landscape context for species richness of bees and wasps in fragmented orchard meadows. Conservation Biology 17:1036–1044.

ter Braak, C. J. F. 1986. Canonical correspondence analysis: a new eigenvector technique for multivariate direct gradient analysis. Ecology 67:1167–1179.

ter Braak, C. J. F. 1995. Ordination. Pages 91-173 in Jongman, C. J. F. ter Braak, and v. Tongeren, editors. Data Analysis in Community and Landscape Analysis. Cambridge University Press, Cambridge, UK.

ter Braak, C. J. F. and Smilauer, P. 1997. Canoco for Windows Version 4.02. Centre for Biometry, Wageningen, The Netherlands.

ter Braak, C. J. F. and Verdonschot, P. F. M. 1995. Canonical correspondence analysis and related multivariate methods in aquatic ecology. Aquatic Sciences 57:1015–1621.

Wagner, H. H. and Fortin, M.-J. 2005. Spatial analysis of landscapes: concepts and statistics. Ecology 86:1975–1987.

TABLE A1. Patch, site and landscape variables used in the ordination and regression analyses. * = units are an arbitrary, 5-step scale devised by Leathwick et al. (2003), † = variable removed due to multicollinearity (Appendix B).

TABLE A2. Exploratory analysis using DISTLM v.2 to identify the most appropriate dissimilarity matrix to describe beetle community composition. DISTLM v.2 conducts a non-linear multivariate regression analysis of user-defined distance matrices using random permutations of the observations (McArdle and Anderson 2001, Anderson 2003). For each distance matrix, the regression tested for the effects of fragment area, distance to edge and their interaction on beetle species composition. Partitioning of the explained variation is shown only for the Euclidean distance measure, which explained the greatest amount of variance in species composition with respect to the hypothesis tested. Codes for variables are as in Table A1.

TABLE A3. (a) Results of the forward selection procedure to determine which of the 64 environmental variables in the pRDA analysis explained significant variation in beetle species composition among sites. (b) Interset correlations between the 26 environmental variables explaining significant variation in beetle species composition between sites and each of the three canonical pRDA axes. Variables are ordered from most to least significant in the forward selection procedure (all with P < 0.05). λ represents the additional variance explained by environmental variables as they are sequentially added into the model. Correlations in bold are significant at a Bonferroni corrected P value of 0.002. Codes for environmental variables are as in Table A1.

FIG. A1. Map of the Hope River Forest Fragmentation Project, South Island, New Zealand. For clarity, small fragments less than 3 ha in area are highlighted with open circles. The continuous forest area extends north and west of the Hope River region to the coastline.

FIG. A2. (a) Flight intercept trap (FIT) used in the Hope River Forest Fragmentation Porject. The FITs were dug into the ground in the forest (b) and matrix (c) sites to act as a combined FIT-pitfall trap.

FIG. A3. Results of pRDA ordination on beetle species composition at 233 sites with respect to 91 patch, site and landscape variables (for full list of variables and species see Table A1 and Appendix C, respectively).

FIG. A4. The slope of the relationship between species diversity and habitat area depends sensitively on the distance from edge at which sampling is conducted. Mean (± 1 SE) slope values for diversity-area curves change non-linearly depending on how far from the forest edge beetles are sampled, when the diversity is represented by both (a) Fisher’s α and (b) Shannon diversity. Effect strengths are measured relative to zero (dashed line); values of zero indicate no effect and deviations from zero (either positive or negative) indicate a strengthening of the effect. Fitted curves are polynomial regressions calculated from data weighted by sample size. Formulae for the fitted lines are as follows, with S = strength of area effect E = distance to edge: (a) S = 0.0668 + 0.0482 × E + 0.0052 × E2; and (b) S = 3.913 + 2.592 × E + 0.274 × E2.