Rebecca E. Campbell, Jon S. Harding, Robert M. Ewers, Stephen Thorpe, and Raphael K. Didham. 2011. Production land use alters edge response functions in remnant forest invertebrate communities. Ecological Applications 21:3147–3161.

Appendix A. A map and detailed attributes of the study sites, information on the sampling design, invertebrate trapping methods, a list of the relative abundances of taxonomic classes and orders captured, a list of environmental and landscape variables and how they were quantified, statistical methodology, formulas and best-fit models for the edge response functions in the manuscript, and output results from the ordination procedures.

Supplementary Methods

Study area and sampling design

The study sites were located on Banks Peninsula (43°43’ – 43°48’S, 172°41’ – 172°54’E) on the east coast of the South Island of New Zealand (Fig. A1). Native forest patches were carefully and objectively matched on the basis of pre-defined criteria, with the aim of reducing the potential confounding effects of variation in fragment location, shape attributes, altitude, aspect, and orientation of the edge (Table A1, Fig. A2).

Invertebrate sampling

The flight intercept traps (FITs) were based on a design by Basset (1988) and consisted of a 40 × 40 cm plastic rain cover, two 60 × 23 cm clear polycarbonate sheets interlocked at right-angles to each other, above a 23 cm diameter funnel and collecting jar (Fig. A3). The funnel and collecting jar were dug into the ground so that the rim of the funnel was flush with the surface of the soil and leaf litter. This allowed the collection of both ground-running and low-flying invertebrates. Care was taken not to disturb the leaf litter or soil immediately surrounding the traps when digging the hole in the ground. The hole was then lined with a metal sleeve so that the trap could be removed and the samples retrieved multiple times without disturbing the surrounding soil further. When the traps were in place, 100 ml of ethylene glycol and 100 ml of absolute ethanol were placed in the collecting jar of each trap to preserve captured invertebrates. Traps were placed a minimum of 15 trap-heights apart from each other (i.e. > 10 m), to minimise potential trap interference for traps at consecutive sampling distance close to the forest edge. This minimum spacing was achieved by offsetting traps laterally while maintaining the required perpendicular distance to the edge.

At the time of collection, samples were sieved through a 125 µm sieve, placed in plastic whirlpak® bags containing absolute ethanol as a preservative and returned to the laboratory for further sorting.

Measured predictors of beetle community composition

Study design variables

For habitat contrast variables (PATCH, MATRIX, NGPATCH, NGMATRIX, NPPATCH, NPMATRIX), the FIT at the 0 m sampling location was considered to be part of the PATCH variable.

Patch area (AREA, ha) was determined from a Geographic Information System (GIS) analysis of the New Zealand Land Cover Database (LCDBv2, Thompson et al. 2003), but we also ground-truthed all sites to verify spatial patch attributes. As a result, we manually re-calculated native forest patch area for Okuti SE (Reynolds1 and Reynolds 2) due to discrepancies between the LCDBv2 and measurements on the ground (e.g. Reynolds 2 was separated by a road which was wide enough to separate the canopies of two native forest patches by greater than 3 m, creating a physical edge). For native patches adjacent to pine plantations at Kinloch Rd in Okuti NW and Middle Rd in Pigeon Bay we could find no record of native forest cover on the LCDBv2, so we calculated the area of the patches manually.

Local-scale biotic and abiotic correlates

Microclimate - Average temperature (TEMP, °C) and relative humidity (RH, %) and incident light were measured using hand held digital meters. Variables were recorded at approximately 50 cm above the ground surface at each of the 72 sampling sites over the period from 28 January to 14 February 2006. For any given patch-to-matrix gradient, microclimate was measured at all nine edge distances on the same day, with readings taken multiple times while walking continuously back and forth along the gradient between 10:00 and 15:00 h. As we were not primarily interested in the absolute microclimate values themselves, we converted all measurements to relative differentials by subtracting the observed value at each sampling distance from the average value recorded in the nearest (within < 1 km) open grassland habitat over the same time interval.

Vegetation - Vegetation structure was recorded within one 5 × 5 m plot centred on each of the 72 FIT traps. Plant species composition and percent cover for each species were estimated visually in three height categories (ground: < 0.5 m; shrub: 0.5–1.5 m; and sub-canopy: 1.5–8 m). Percent cover scores within tiers were converted into a single importance value for each species by weighting individual tier scores by their corresponding tier depths, using the formula:

Subsequently, a Detrended Correspondence Analysis (DCA), with down-weighting of rare species, was used to convert the single importance values for each plant species to an overall measure of variation in plant community composition. The first three DCA axes (PLANT1, PLANT2 AND PLANT3) were used as vegetation correlates of variation in beetle community composition. The DCA was conducted in CANOCO (Version 4.02) (ter Braak and Smilauer 1997).

Leaf litter complexity - We randomly placed five 10 × 10 cm leaf litter quadrats within each of the 5 × 5 m vegetation plots. Within each of the 10 × 10 cm quadrats we collected all the leaf litter down to the compact soil layer, placed it in an individual plastic bag, and returned it to the laboratory where it was frozen to preserve the samples until they could be processed. We oven-dried the litter for 72 hours at 50°C and recorded total litter dry mass (LITTERMASS, g/m²). We then sieved the litter into six different size classes: < 790 µm, < 1.4 mm, < 10 mm, < 12.7 mm, < 20 mm and ≥ 20 mm (cylindrical mesh dimension in the latter two cases, and square mesh dimension in the others). We individually weighed each size class fraction and calculated leaf litter complexity (LITTERDIV) using the Shannon-Weiner diversity index:

where p_i is the proportion of litter mass in each size class category relative to the total number of categories. Lastly, we obtained a measure of litter composition from axis 1 scores of a DCA ordination of the relative mass of litter in the six size classes (LITTER1).

Invertebrate ordinal abundance - A DCA ordination, with down-weighting of rare taxa, was conducted on the relative abundances of classes and orders of invertebrates, excluding beetles, in order to determine whether the responses of other invertebrate taxa might be correlated with beetle species responses to land use change. Axes 1, 2 and 3 of the DCA ordination (INVERT1, INVERT2, INVERT3) were included as environmental correlates in the species-level analysis of beetle community composition.

Patch-scale spatial correlates

We obtained the absolute perimeter length of patches (PERIM, m) from a GIS analysis of the LCDBv2, but we also ground-truthed all sites as described above for AREA, and recalculated PERIM accordingly.

To derive an area-standardised measure of patch shape complexity, we calculated the shape index (SI) of Patton (1975):

where PERIM is in metres and AREA is in hectares.

Landscape-scale land cover correlates

Using the standard buffering and clipping techniques available in ArcGIS, we created circular buffers (‘landscapes’) of 243 m, 729 m and 2178 m radius for the LCDBv2 land cover data around each of the 72 trap locations. The spatial scale of these ‘landscapes’ was selected as a logical extension of the log3-scale edge sampling distances, and encompassed the range of spatial scales over which invertebrates are thought to respond to landscape structure (e.g., Holland et al. 2004). Within each landscape, five variables were calculated representing variation in landscape composition: the number of native forest patches (NATFRAG243, NATFRAG729 and NATFRAG2178, for each of the three spatial scales, respectively), percent cover of native woody vegetation (NAT243, NAT729, NAT2178), percent cover of pine plantation (PINE243, PINE729, PINE2178), percent cover of grassland pasture (GRASS243, GRASS729, GRASS2178), and land cover diversity, calculated using the Shannon-Wiener index (DIV243, DIV729, DIV2178). Native woody vegetation cover included ‘grey scrub’ (mostly Coprosma spp.), manuka/kanuka scrub, and native forest (LCDBv2).

Potential confounding variables

Although sample effort was approximately consistent between sites, and invertebrate capture rates were standardised for number of trap days sampling effort, there is still the potential for total sample abundance to have a confounding effect on comparisons of community composition if abundance is highly variable across sites. Consequently, sample abundance (BEETLEN) was included as a potential confounding variable in multivariate analyses.

Although every effort was made to control variation in patch attributes during site selection, minor variation in the orientation of the sampled edge (ORIENT, °), the aspect of the native forest patch (ASPECT, °) and the average elevation of the patch (ALT, m), might represent potential confounding factors in the analysis. Edge orientation was recorded as positive or negative deviations from north (0°). Aspect and elevation were determined from the New Zealand Digital Elevation Model (DEM) using ArcGIS.

Land-use treatments were paired within catchment locations as part of the spatial blocking in the study design, but as we were not interested in the absolute differences between catchments in this study, we treated them as potential confounding variables in the analysis (PIGEON, KAITUNA, OKUTI), in order to remove any spatial autocorrelation in beetle community composition attributable to variation among catchments.

To deal with any residual spatial autocorrelation among sampling locations (e.g. within catchments), we used a cubic trend surface approach with nine linear, quadratic and cubic combinations of latitude and longitude coordinates (LAT, LONG, LATLONG, LAT2, LONG2, LAT2LONG, LAT LONG2, LAT3, LONG3) from the NZ Metric Grid (NZMG). The NZMG values were recoded to minimum northings and eastings within the study area, and truncated to a 1000 m resolution so as to avoid removing the intrinsic fine-scale spatial autocorrelation within edge gradients that we were interested in measuring.

SUPPLEMENTARY ANALYSES

Standardising measures of abundance and species richness

Invertebrate trapping effort varied from 13–22 days sampling across the 72 sites. Standardisation for differential trapping effort involved conversion of sample abundance to a flux density estimate, calculated as the number of individuals per m² of trap surface per trap-day. The (single-sided) surface area of the FIT trap was 0.276 m².

Sample-based species accumulation curves were calculated using EstimateS 7.5.0 (Colwell 2005), with estimated species richness and the frequency of singletons and doubletons (species only represented by one or two individuals, respectively, in a given random draw) expressed against sample abundance re-scaled by number of individuals. Accumulation curves were calculated with 100 replicate random draws (without replacement). The Chao1 species richness estimator was calculated to compare predicted, asymptotic species richness (at high sample abundance) across production land use treatments. For comparisons of land-use type (native forest versus pasture versus pine plantation), species accumulation curves were calculated without the 0 m edge distance (which cannot be assigned unequivocally to one land use type). For comparisons of the effect of production land use on estimated species richness within native forest patches (native patches adjacent to pasture versus native patches adjacent to pine), species accumulation curves were calculated with the inclusion of the 0 m edge distance (because in this case species turnover at 0 m is a component of the matrix influence being tested).

Species rarefaction was conducted in Biodiversity Professional software (BDPro v3.2) (McAleece 1997) and involved an individual-based randomisation procedure in which the number of species at each site was estimated after first standardising to the lowest common sample abundance across sites. We chose a conservative rarefaction cut-off (knot) of N = 41, with sites having total sample abundance smaller than this cut-off being dropped from the analysis. Samples that were excluded were 7 (out of a total of 16) samples from pasture habitat, 2 (out of 16) from pine habitat, 2 (out of 20) from native fragments adjacent to pine and 1 (out of 20) from native fragments adjacent to pasture. A lower cut-off value would have included more sites, but would have rapidly weakened the value of the comparison, as the low abundance of the random draw constrains the range of possible variation in species richness.

Constrained, direct-gradient analysis of variation in beetle community composition

Canonical Correspondence Analysis (CCA) is a direct-gradient multivariate ordination technique that relates community composition to known variation in the environment (ter Braak 1986). It assumes the existence of a set of underlying environmental gradients that all species respond to, and provides an estimate of species-environment relationships (ter Braak 1986).

In an initial partial CCA analysis (pCCA), we accounted for variance due to the experimentally-controlled blocking effect of river catchments by entering PIGEON, KAITUNA, and OKUTI as covariables, and then tested whether any of the remaining 13 potential confounding variables explained significant variance in the data. We used the forward selection procedure (with 9999 permutations) to rank potential confounding variables in order of their importance in explaining variation in species composition (ter Braak 1986). In the forward selection procedure, the first variable selected is the variable with the highest marginal eigenvalue (fit with only one variable in the model), and subsequent variables are selected with the highest successive conditional eigenvalues (additional fit when further variables are added) (Jongman et al. 1995, Leps and Smilauer 2003). Only two further confounding variables (ALT and ASPECT) explained significant variation in beetle species composition, and the remaining 11 variables were discarded. Therefore, in the final pCCA there were five covariables entered in the model (PIGEON, KAITUNA, OKUTI, ALT, ASPECT) prior to testing the effects of the 15 treatment variables and 29 local-, patch- and landscape-scale environmental correlates (again, using the forward selection procedure, as described above).

The sample scores used for the pCCA biplot graphs were linear combinations of environmental variables (rather than weighted average scores) in order to better visualise similarities and differences in community composition with respect to the significant environmental variables (Leps & Smilauer 2003). Significant environmental variables were plotted as vectors on the pCCA plot. Site centroids and the centroids for several beetle species which were significantly correlated with the pCCA axes were plotted to show (in an approximate way) their distribution in relation to each environmental variable (ter Braak 1986).

Intraset correlations were calculated to show the relationship between the environmental variables and site ordering along the pCCA axes (Jongman et al. 1995) and correlations between individual beetle species and pCCA axes 1 and 2 scores were used to determine which species were significantly correlated with the major axes among sites.

Quantifying edge response functions of invertebrates

For each of four major categories of invertebrate response variables (capture rates of invertebrate classes and orders, beetle species richness, variation in beetle community composition, and capture rates of individual beetle species), continuous edge response functions were calculated across patch-to-matrix gradients for each of the two contrasting adjacent production land uses, using the statistical approach of Ewers and Didham (2006). Using a form of the general logistic model we determined the best-fit edge model out of five models of increasing complexity (Ewers and Didham 2006):

(1) the null hypothesis of no discernable edge effect, calculated as the mean of the response variable η:

where ε is an error term;

(2) a simple linear gradient of the form:

where β₀ and β₁ are constants and D is the distance to edge;

(3) a power model:

(4) a logistic model that describes a sigmoidal change in community composition across an edge, with an asymptote in both the patch and matrix habitats:

which also includes an additional constant, β₃;

(5) a unimodal model based on the logistic model, but with one extra constant (β₄) and a D² term to describe a unimodal change in community composition at a particular distance from an edge:

We fitted these five models to patch-matrix gradients, treating the four patch gradients within each of the two production land-use treatments (pasture matrix versus pine matrix) as replicates. In each case, we assessed model significance and calculated the Akaike Information Criterion (AIC) value for each model. We selected the best model as being the one with the lowest AIC value, or in the case of multiple models within two AIC units of each other, we selected the simplest model (with the fewest parameters). Model fitting was conducted in R version 2.5.1 (R Development Core Team 2004).

When determining the edge response functions of multiple invertebrate classes or orders, and multiple individual beetle species, the total abundance and frequency of occurrence of taxa across the 72 sites had a major bearing on model fit, and the interpretation of edge responses. Therefore, we conducted a post-hoc power analysis in R to determine the minimum sample abundance necessary to detect a large effect size of 0.35 (sensu Cohen 1988) for the most complex regression model (the unimodal model with 5 parameters), using conventional values of power = 0.8 and α = 0.05. To meet these criteria, total sample abundance would have to be N > 42.5, therefore we only tested our edge response functions for taxa with 43 or more individuals.

SUPPLEMENTARY RESULTS

No variation in beetle species richness with land use

In order to assess sampling adequacy, a sample-based rarefaction analysis was used to plot estimated species richness against sampling intensity, but the resulting species accumulation curve for all sites combined did not asymptote (Fig. A4), indicating that our sampling regime did not collect all species present in the regional species pool. However, we conservatively estimate that we sampled between 65.1 % to 82.8 % (calculated from Chao 1 index values in EstimateS) of the total number of species likely to be captured using this trapping method across native, pine and pasture land uses on Banks Peninsula. The number of species represented by only one or two individuals (singletons and doubletons) stabilised with increasing sample effort, but did not decrease even at the maximum sampling effort employed (Fig. A4).

Beetle capture rates varied substantially between land use types, but this trend did not arise as a sampling artefact, as beetle abundance per trap was not significantly correlated with sampling effort (number of trap days) for this collection period (R = 0.1997, N = 72, P = 0.093).

TABLE A1. Patch and landscape attributes of the study sites, selected based on objective pre-defined criteria (see details in text) so that average characteristics were as similar as possible between the two categories of land-use contrast: native forest patches adjacent to exotic pasture and native forest patches adjacent to pine plantations, paired within catchment locations, Banks Peninsula, South Island, New Zealand.

* data from Eikaas et al. (2005)
^§ Within a radius of 2187 m around the native patch, calculated using ArcGIS on LCDB2 land cover
^†Calculated using ArcGIS on LCDB2 land cover
^‡Calculated using ArcGIS and the New Zealand Digital elevation model (DEM).

TABLE A2. Description of study design variables, and local-scale, patch-scale, and landscape-scale environmental variables tested as predictors of beetle community composition in the forward selection procedure of a canonical correspondence analysis (CCA) ordination, as well as variables considered to have potential confounding effects on the interpretation of patterns of variation in beetle community composition. NZMG = New Zealand Metric Grid mapping coordinate system.

TABLE A3. Total abundance of terrestrial invertebrate taxa captured using flight interception traps at nine distances from the edge of each of four replicate native forest patches adjacent to a low-contrast pine plantation matrix, and four replicate native forest patches adjacent to a high-contrast pasture matrix on Banks Peninsula, Canterbury, New Zealand. Abundance totals are tabulated for the five distances int the native forest patches (including the 0 m edge trap), and four distances into the respective matrix habitats, with the total trap days sampling effort indicated. ‘Other’ represents larvae and pupae, which were not assigned to taxonomic group.

TABLE A4. A comparison of Akaikie Information Criterion (AIC) values used to identify the best-fit continuous edge response function (indicated in bold) from a family of five competing models of increasing complexity. Edge response functions were calculated separately for (a) the four replicate native forest patches adjacent to pine matrix, and (b) the four replicate native forest patches adjacent to pasture matrix. Response variables tested were the capture rates of Araneae (see Fig. 1a), Coleoptera (see Fig. 1b), Diplopoda (see Fig. 1c), Hemiptera (see Fig. 1d), the absolute species richness of all Coleoptera (see Fig. 2a), the richness of Coleoptera standardised for differing samples abundance (see Fig. 2b), variation in Coleoptera community composition (see Fig. 4), the capture rates of sp. 2276 Melanophthalma discoidea and sp. 2278 Corticariinae sp. 6 (both Latridiidae: Corticariinae), and the capture rates of sp. 2317 Acrotrichis sp. and sp. 2318 Notoptenidium sp. (both Ptiliidae) (see Fig. 5a–d). ‘Community composition’ refers to the axis 1 scores of a partial Canonical Correspondence Analysis ordination, as depicted in Fig. 3. NA indicates the model fitting procedure did not converge. For models with AIC values within two units of each other, the model with the fewest parameters was selected as the best model. A dash indicates that there were no specimens collected in this land use.

TABLE A5. Results of a partial Canonical Correspondence Analysis of variation in beetle species composition in relation to measured environmental heterogeneity and land use change. Analyses were conducted using a two-stage partialling procedure, first taking into account the spatial ‘pairing’ of treatment replicates imposed within catchments (PIGEON, KAITUNA, OKUTI), and second taking into account the confounding effects of elevation (ALT) and aspect (ASPECT). Eigenvalues are a measure of the variance in community composition explained by each of the first four orthogonal axes of total variation (i.e. total inertia) in the data.

FIG. A1. Location of the native forest fragments sampled on Banks Peninsula, South Island, New Zealand. The map shows the dominant surrounding land cover types from the New Zealand Land Cover Database (LCDB v.2) in the Kaituna, Pigeon Bay and Okuti catchments.

FIG. A2. Following careful selection of native forest patches based on pre-defined objective criteria (see text), patch characteristics did not vary significantly between the two land-use contrast categories: (a) patch area, (b) patch elevation, and (c) shape index. See Table A1 for statistical significance.s

FIG. A3. The flight interception trap design, modified from Basset (1988), illustrating how the collecting funnel and jar were dug into the ground, and the hole was lined with a tin sleeve so that the top of the funnel opening was flush with the ground surface, allowing the capture of both surface active and flying invertebrates. Dimensions are as described in the text (figure reproduced with permission from Chris Edkins, Department of Conservation).

FIG. A4. Species accumulation curves were calculated (a) for all 72 sampling sites combined, (b) for the three differing land use types separately, and (c) for native forest fragments adjacent to differing matrix land use types. Mean (± 1SE) estimated species richness was calculated using EstimateS v.7.5 with 100 random draws, and samples drawn without replacement. ‘Singletons’ indicates the number of species with a single individual from a random sample of the indicated sample size; similarly, ‘doubletons’ indicates the number of species with two individuals. Calculations in (b) exclude the zero metre edge distances which are difficult to objectively assign to land use type, whereas calculations in (c) do include the edge (0 m) distance, as this is a representation of how species mix from the directly adjacent land use.

LITERATURE CITED

Cohen, J. 1988. Statistical Power Analysis for the Behavioral Sciences. Lawrence Erlbaum Associates, Hillsdale, NJ.

Colwell, R. K. 2005. EstimateS. Statistical estimation of species richness and shared species from samples. Version 7.5. URL: [purl.oclc.org/estimates].

Eikaas, H. S., A. R. McIntosh and A. D. Kliskey. 2005. Catchment- and site-scale influences of forest cover and longitudinal forest position on the distribution of a diadromous fish. Freshwater Biology 50:527–538.

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

Ewers, R. M., and R. K. Didham. 2006. Continuous response functions for quantifying the strength of edge effects. Journal of Applied Ecology 43:527–536.

Holland, J. D., D. G. Bert, and L. Fahrig. 2004. Determining the spatial scale of species' response to habitat. BioScience 54:227–233.

Jongman, R. H. G., C. J. F. ter Braak, and O. F. R. van Tongeren, editors. Data Analysis in Community and Landscape Analysis. Cambridge University Press, Cambridge.

Leps, J. and P. Smilauer. 2003. Multivariate Analysis of Ecological Data using CANOCO. Cambridge University Press, Cambridge.

McAleece, N. 1997. BioDiversity Pro. The Natural History Museum and Scottish Association for Marine Science, London, UK.

R Development Core Team. 2004. R Foundation for Statistical Computing, Vienna.

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., and P. Smilauer. 1997. Canoco for Windows Version 4.02. Centre for Biometry, Wageningen, The Netherlands.

Thompson, S., I. Grüner, and N. Gapare. 2003. New Zealand Land Cover Database Version 2: Illustrated guide to target classes. New Zealand Ministry for the Environment., Wellington, New Zealand.