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 lowlevel intercept and pitfall trap. This trap type was chosen because it is suitable for capturing both grounddwelling 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 fragmentlevel 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_{2} 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^{2}) 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 ShannonWeiner 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 7d 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 5min 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 5min 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 5min 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 5min 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 5min light reading was averaged with the three earlier and three later 5min values to give a 30min slidingwindow smoothing function. Variability in site light environment (PFDVAR) was measured as the CV of the 5min 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 5min 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 5min readings were missing within any 4hr 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.
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), ShannonWeiner 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 (SteffanDewenter 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 fourlevel 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 speciesenvironment 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 GarciaBoyero 1995, AbenspergTraun and Smith 1999, Kotze and Samways 1999, 2001, Davies et al. 2001, Barbosa and Marquet 2002, SteffanDewenter 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, AbenspergTraun 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 wholecommunity 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_{10}(log_{10}(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 oneway 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 nonpreferred 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_{10}(log_{10}(A) + 2) + 0.25 × (log_{10}(log_{10}(A) + 2)^{2}; and (b) S = 0.067 + 0.048 × log_{10}(log_{10}(A) + 2) + 0.005 × (log_{10}(log_{10}(A) + 2)^{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_{2}(│E│)) + 0.002 × (log_{2}(│E│))^{2} ; and (b) S = 0.0736 + 0.1274 × (log_{2}(│E│)) + 0.0117 × (log_{2}(│E│))^{2} .
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
AbenspergTraun, M., and G. T. Smith. 1999. How small is too small for small animals? Four terrestrial arthropod species in differentsized remnant woodlands in agricultural Western Australia. Biodiversity and Conservation 8:709–726.
AbenspergTraun, M., G. T. Smith, G. W. Arnold, and D. E. Steven. 1996. The effects of habitat fragmentation and livestockgrazing 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 distancebased 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. GarciaBoyero. 1995. The effects of forest fragmentation on butterfly communities in central Spain. Journal of Biogeography 22:129–140.
Borcard, D., P. Legendre, C. AvoisJacquet, 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 betweenpatch 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 woodlandpasture boundaries in the central wheatbelt 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 distancebased redundancy analysis. Ecology 82:290–297.
Pajunen, T., Y. Haila, E. Halme, J. Niemelä, and P. Punttila. 1995. Grounddwelling spiders (Arachnidae, Araneae) in fragmented old forests and surrounding managed forests in southern Finland. Ecography 18:62–72.
SteffanDewenter, I. 2002. Landscape context affects trapnesting bees, wasps, and their natural enemies. Ecological Entomology 27:631–637.
SteffanDewenter, 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 91173 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, 5step scale devised by Leathwick et al. (2003), † = variable removed due to multicollinearity (Appendix B).


Code 
Explanation 
Units 
Patch attributes 




DISTCODE 
log_{2} Distance to edge 
± m 


ABSDIST 
absolute value of DISTCODE 
m 


LOGAREA 
Log_{10} Area 
ha 

† 
LOGSHAPE 
log_{10} Shape Index 



AREAEDGE 
Areaedge interaction (LOGAREA × DISTCODE) 



AREAABSD 
Areaabsolute edge distance interaction (LOGAREA × ABSDIST) 



FORCONT 
Deep forest dummy variable 
binary 


MATCONT 
Deep matrix dummy variable 
binary 





Site attributes 




BIOMASS 
Total basal area of standing live trees 
m^{2}/ha 


CVTREE 
Coefficient of variation of tree dbh 
cm 


TREEDIV 
Tree species richness 



VEGPCO1 
Axis 1 from PCoA of tree composition 



VEGPCO2 
Axis 2 from PCoA of tree composition 



DEADST 
Total volume of standing dead wood 
m^{3}/ha 


CVDEADST 
Coefficient of variation of volume of standing dead wood 
% 


WOOD 
Total volume of coarse woody debris 
m^{3}/ha 


CVWOOD 
Coefficient of variation of volume of coarse woody debris 
% 


CANOPY 
Canopy height 
m 


LEAFLIT 
Mean dry weight mass of litter 
kg/ha 


CVLEAF 
Coefficient of variation of five litter samples 
% 


QDIV 
Shannon diversity index for quadrat ground cover categories 



QUADPCO1 
Axis 1 score from PCoA of quadrat ground cover 



QUADPCO2 
Axis 2 score from PCoA of quadrat ground cover 



QUADPCO3 
Axis 3 score from PCoA of quadrat ground cover 


Microclimatic variables 




DELTAT 
Average daily temperature differential between site and 
°C 


TVAR 
Average daily temperature range 
°C 


DELTARH 
Average daily relative humidity differential between site and openarea reference point 
% 


RHVAR 
Average daily coefficient of variation of relative humidity 
% 


DELTAPFD 
Average daily incident light 
% 


PFDVAR 
Average daily coefficient of variation of incident light 
% 





Landscape attributes 




NINADIST 
Distance to continuous forest 
m 


256NUM 
Number of forest fragments within radius of 256 m 



512NUM 
Number of forest fragments within radius of 512 m 



1024NUM 
Number of forest fragments within radius of 1024 m 



2048NUM 
Number of forest fragments within radius of 2048 m 


† 
4096NUM 
Number of forest fragments within radius of 4096 m 



8192NUM 
Number of forest fragments within radius of 8192 m 



256FOR 
Forest cover in landscape with radius of 256 m 
% 

† 
512FOR 
Forest cover in landscape with radius of 512 m 
% 


1024FOR 
Forest cover in landscape with radius of 1024 m 
% 

† 
2048FOR 
Forest cover in landscape with radius of 2048 m 
% 


4096FOR 
Forest cover in landscape with radius of 4096 m 
% 


8192FOR 
Forest cover in landscape with radius of 8192 m 
% 


CON256 
Forest connectivity index for landscape with radius of 256 m 



CON512 
Forest connectivity index for landscape with radius of 512 m 



CON1024 
Forest connectivity index for landscape with radius of 1024 m 



CON2048 
Forest connectivity index for landscape with radius of 2048 m 


† 
CON4096 
Forest connectivity index for landscape with radius of 4096 m 


† 
CON8192 
Forest connectivity index for landscape with radius of 8192 m 



DIV256 
Landscape diversity for landscape with radius of 256 m 



DIV512 
Landscape diversity for landscape with radius of 512 m 



DIV1024 
Landscape diversity for landscape with radius of 1024 m 



DIV2048 
Landscape diversity for landscape with radius of 2048 m 



DIV4096 
Landscape diversity for landscape with radius of 4096 m 


† 
DIV8192 
Landscape diversity for landscape with radius of 8192 m 


Soil composition 





ACIDP 
Acid soluble phosphorus content of soil 
% 


CALCIUM 
CALCIUM content of soil 
% 

† 
PSIZE 
Soil particle size 
* 

† 
AGE 
Soil age 
* 

† 
INDUR 
Soil parent material hardness 
* 

Macroclimate 





JUNES 
Winter solar radiation 
MJ.m^{2}.d^{1} 


MAS 
Annual solar radiation 
MJ.m^{2}.d^{1} 

† 
MAT 
Mean annual temperature 
ºC 


TMIN 
Winter minimum temperature 
ºC 


R2PET 
Water balance ratio 



DEFICIT 
Annual water deficit 
mm 

† 
VPD 
October vapour pressure deficit 
kPa 

Land Environments of New Zealand categories 




P 
LENZ level 1 dummy variable 
binary 


K 
LENZ level 1 dummy variable 
binary 

† 
P1 
LENZ level 2 dummy variable 
binary 

† 
K1 
LENZ level 2 dummy variable 
binary 


E1 
LENZ level 2 dummy variable 
binary 


K11 
LENZ level 3 dummy variable 
binary 

† 
P12 
LENZ level 3 dummy variable 
binary 


E14 
LENZ level 3 dummy variable 
binary 


K11D 
LENZ level 4 dummy variable 
binary 


E14C 
LENZ level 4 dummy variable 
binary 


E42A 
LENZ level 4 dummy variable 
binary 


K11B 
LENZ level 4 dummy variable 
binary 





Confounding variables 




ALTITUDE 
Altitude 
m 


SLOPE 
Slope 
º 


ASPECT 
Aspect 
º 


LAT 
Latitude 
1000 m 


LONG 
Longitude 
1000 m 


LATLONG 
Latitudelongitude spatial autocorrelation variables 



LAT2 
Latitudelongitude spatial autocorrelation variables 



LONG2 
Latitudelongitude spatial autocorrelation variables 



LAT2LONG 
Latitudelongitude spatial autocorrelation variables 



LATLONG2 
Latitudelongitude spatial autocorrelation variables 



LAT3 
Latitudelongitude spatial autocorrelation variables 



LONG3 
Latitudelongitude spatial autocorrelation variables 

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 nonlinear multivariate regression analysis of userdefined 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.
Distance measure 
df 
SS 
MS 
F 
P 
r^{2} 

BrayCurtis 








Regression 
3 
20.68 
6.89 
25.33 
< 0.001 
0.25 

Residual 
229 
62.31 
0.27 











Chi square 








Regression 
3 
277.63 
92.54 
4.68 
0.003 
0.06 

Residual 
229 
4531.30 
19.79 











Euclidean 








Regression 
3 
3225.18 
1075.06 
26.04 
< 0.001 
0.25 

LOGAREA 
1 
999.17 
999.17 
24.20 
< 0.001 
0.08 

DISTCODE 
1 
1796.21 
1796.21 
43.50 
< 0.001 
0.14 

AREADIST 
1 
429.80 
429.80 
10.41 
0.001 
0.03 

Residual 
229 
9455.96 
41.29 



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.

(a) Forward selection 

(b) Interset correlations 

Variable 
λ 
P 
F 

Axis 1 
Axis 2 
Axis 3 








CANOPY 
0.1 
0.001 
36.06 

0.8491 
0.1098 
0.0716 
QUADPCO1 
0.02 
0.001 
5.9 

0.8344 
0.1051 
0.0752 
MATCONT 
0.01 
0.001 
5.02 

0.1827 
0.0164 
0.5621 
AREADIST 
0.01 
0.001 
4.83 

0.6738 
0.2843 
0.1149 
CON512 
0.01 
0.001 
4.31 

0.0159 
0.3852 
0.2041 
CON1024 
0.01 
0.001 
4.26 

0.2559 
0.2784 
0.021 
LOGAREA 
0.02 
0.001 
4.19 

0.1204 
0.281 
0.0297 
CON2048 
0.01 
0.001 
4.17 

0.0825 
0.168 
0.1892 
DISTCODE 
0.01 
0.001 
4.08 

0.8072 
0.0411 
0.2544 
ABSDIST 
0.01 
0.001 
3.1 

0.0751 
0.2259 
0.2591 
CON256 
0 
0.001 
2.63 

0.1827 
0.0724 
0.3565 
NINADIST 
0.01 
0.001 
2.29 

0.045 
0.0517 
0.4582 
QUADPCO3 
0.01 
0.001 
2.16 

0.1209 
0.3241 
0.0628 
CALCIUM 
0 
0.001 
2.1 

0.2492 
0.0383 
0.0054 
FORCONT 
0 
0.001 
2.03 

0.0177 
0.0349 
0.3605 
AREAABSD 
0.01 
0.001 
2.01 

0.0559 
0.3319 
0.139 
DELTARH 
0.01 
0.001 
1.97 

0.6676 
0.1637 
0.0539 
TVAR 
0.01 
0.001 
1.75 

0.7193 
0.0607 
0.1004 
TREEDIV 
0.01 
0.001 
1.65 

0.177 
0.1336 
0.0248 
LEAFLIT 
0 
0.002 
1.57 

0.101 
0.0294 
0.0186 
E14 
0 
0.002 
1.52 

0.3729 
0.0432 
0.1748 
MAS 
0 
0.006 
1.47 

0.0409 
0.2019 
0.0623 
8192FOR 
0 
0.008 
1.55 

0.1554 
0.0816 
0.178 
8192NUM 
0.01 
0.01 
1.5 

0.0407 
0.0487 
0.2017 
CVTREE 
0 
0.035 
1.36 

0.7133 
0.1586 
0.0237 
QDIV 
0 
0.049 
1.28 

0.0151 
0.1232 
0.0422 
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 FITpitfall 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 diversityarea curves change nonlinearly 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 × E^{2}; and (b) S = 3.913 + 2.592 × E + 0.274 × E^{2}. 