Species-specific characteristics such as body size or dispersal ability are likely to be crucial in determining which aspects of habitat heterogeneity e.g. availability, patchiness, structural diversity are relevant and at which spatial scales individual perceive landscape structure, or are affected by habitat fragmentation (Roland and Taylor, 1997; Steffan-Dewenter et al., 2002; Chust et al., 2003; Dauber et al., 2005; Wiens and Milne, 1989). Furthermore, the response to habitat heterogeneity may also be scale-specific (Vanbergen et al., 2007). For example, carabid beetle assemblages have been shown to be positively influenced by habitat heterogeneity at micro-scales (0.25 m2) and mesoscales (500-1000 m2) but not at macroscales (10 km2) (Brose, 2003; Tews et al., 2004).
The influence of habitat fragmentation on population dynamics have often been confounded with those of habitat composition (i.e. habitat amount and quality; Ewers and Didham, 2006) and habitat loss (Fahrig 2003). Empirical studies have shown that fragmentation might have positive or negative effects on populations, but in order to determine the factors that lead to such consequences, more studies that separate between effects of habitat loss and fragmentation per se are needed (Fahrig, 2003). Spatially-explicit simulation models, in which these two factors can be manipulated independently (e.g., Nabe-Nielsen et al., 2010) can be the key to produce the required evidence.
Neutral landscape models (NLMs), in particular, provide a good modelling framework to study habitat fragmentation. NLMs are grid-based maps in which complex habitat distributions are generated by random, hierarchical, or fractal algorithms. Because the artificial landscapes are generated with analytical algorithms, they are thus "neutral" to the biological and physical processes that shape real landscape patterns. In a NLM cells are identified by habitat type or some other landscape feature; using fractal algorithms to generate NLMs allows simple control over spatial autocorrelation, and complex landscapes can be generated in a systematic way, just by varying the degree of spatial autocorrelation of habitat patches (With, 1997).
Neutral landscape models were introduced to generate spatial patterns in the absence of any structuring process (Hargrove et al., 2002), provide null models for predicting when habitat fragmentation occurs and is expected to affect population dynamics (With, 1997).
Microarthropod communities have often been used as a model system to study effects of habitat fragmentation, and results of several studies indicate their usefulness for this purpose. For instance, microarthropods have been used to investigate how autocorrelation of disturbance events affects time to extinction of populations (Pike et al., 2004), or how habitat loss and fragmentation affects population viability and ecosystem functioning (Astrom and Bengtsson, 2011).
Habitat fragmentation may be the result of natural phenomena such as fire (Wright, 1974) and windfall (Foster, 1980). However, the most important and large-scale cause of habitat fragmentation is the expansion and intensification of human land use (Burgess and Sharpe, 1981).
Among the different human activities that may cause habitat destruction, we focused on agricultural practice, as it is especially relevant for soil organisms. Robertson et al. (1993) showed that cultivation changes the spatial structure of important soil properties, and several field studies have investigated the relationship between soil metal pollution and microarthropod distribution (see e.g. Bengtsson and Rundgren, 1988; Hagvar and Abrahamsen, 1990; Salminen and Haimi, 1999; Gongalsky et al., 2010).
In most studies, heterogeneity in the distribution of toxic metals in soil has been found to be significantly related to spatial changes in the community structure of springtails and mites (Caruso et al., 2009).
In most agricultural practices, both inorganic and organic compounds are applied to protect crops from pest species. Among inorganic chemicals, copper is widely used as fungicide on a number of crops. Due to its long degradation time, sites where copper has been applied for several years can reach very high concentrations in soil, and this causes a permanent loss of habitat quality for soil organisms. In contrast, organic agrochemicals are often easily degradable, but are generally more toxic to microarthropods than metals: their application thus acts as a disturbance event for these communities, killing a smaller or larger fraction of exposed populations, and quickly degrading afterwards.
Different crops and modes of application of agrochemicals can cause different patterns of distribution of the compounds in soil. For instance, a pesticide which is sprayed on homogeneously distributed crops is more likely to have a random distribution in soil, while a treatment in bands (e.g. potato furrows) will have a spatially correlated distribution. Furthermore, another source of disturbance that comes with agricultural land use is the mechanical stress caused by tillage and other uses of machinery. (Maraun et al., 2003) suggested that Collembola are sensitive to mechanical disturbances, and according to several studies reviewed by (Petersen, 2002) collembolan density reported for some cultivated fields was generally low compared with data from natural or semi-natural sites.
Beside chemical and mechanical stressors, populations can be exposed to physical ones, such as periods of drought, especially during summer, to which Collembola have little resistance (Holmstrup, 1997).
In the present study, we implemented a fractal algorithm in a spatially explicit individual-based model, and focused on the effects of different patterns of habitat fragmentation and disturbance events at a microscale on the population dynamics of the collembolan Folsomia candida (Willem). For this purpose we designed simulation experiments in which three different factors interact with individuals in the model: 1) different spatial patterns of copper contamination represent permanent loss of habitat, 2) a realistic implementation of a summer drought period is used to investigate the influence of natural stress on population-level effects of copper and, 3) different levels and numbers of disturbance events represent the other stressors (i.e., pesticide applications, tillage, etc.) to which springtail populations are exposed in the field.
Furthermore, we ran simulations with and without avoidance behaviour, as it has been shown that F. candida avoid copper but do not detect all toxicants (Greenslade and Vaughan, 2003). This allowed us to investigate which effects habitat loss and fragmentation would hypothetically have on populations in cases in which contaminated habitat is not avoided, and allowed us to test whether knowing if a compound is avoided or not would change the estimation of risk posed by the toxicant.
The aim of our study was to explore hypotheses and scenarios to show how different management options could reduce long term risks for soil invertebrates.
More specifically, we tested the following hypotheses: (i) reduction in population abundance caused by a progressive habitat reduction is more than proportional to the habitat loss; (ii) collembolan populations are more affected by habitat loss when the degree of fragmentation of remaining habitat is higher (i.e. patches are spatially uncorrelated or poorly correlated), especially if the percentage of available habitat is low; (iii) if individuals cannot avoid contaminated patches of soil, population-level effects of habitat loss and fragmentation are worse, especially if the percentage of available habitat is low; (iv) reduction of habitat quality and different patterns of disturbance events, caused by intensive agricultural practices, combined with naturally occurring physical stress impact population recovery and may lead to extinction.
Methods
The species used in the simulations is Folsomia candida Willem 1902, which belongs to the order Collembola, suborder Entomobryomorpha, family Isotomidae. This species is used as a standard test organism for toxicity tests: a 28-day reproduction test (International Organization for Standardization, 1999; Organisation for Economic Co-operation and Development, 2009) is included in the refinement options for ecological risk assessment of plant protection products to soil organisms in the EU .
Copper sulphate (CuSO4) was used as a model contaminant to simulate permanent loss of habitat quality: it is proven to cause toxic effects to F. candida survival and reproduction, and to elicit behavioural responses like avoidance (Boiteau et al., 2011). Moreover it is used as fungicide on a variety of crops, is the most widely distributed pollutant among all metals, and therefore relevant for ecological risk assessment.
Full descriptions of both the biology of Folsomia candida and the individual-based model are found in Meli et al. (2013). Therefore, in the following section we give only a brief overview of the model itself, while we focus on the submodels that have been added to the original model in order to test the hypotheses listed in the Introduction.
Individual-based model overview
The purpose of the model is to investigate how populations of F. candida are affected by spatial distribution of toxic contamination in soil, with a special focus on interactions with food availability and local population density (Meli et al., 2013). The model comprises the entities eggs, juvenile and adult female springtails, and grid cells. Springtails are mobile and are characterized by the state variables age (days), position (continuous coordinates), direction of movement, energetic status (days-to-death), cumulative distance (in cm) walked in each hourly time-step, and time (h) spent on contaminated grid cells. Grid cells are characterized by their food level and concentration of toxicant (mg kg-1 soil). The model world is a two-dimensional grid of 100x100 square grid cells representing 1 cm2 of soil. The global environment is characterized by six "seasons", which determine the temperature-dependent life-cycle parameters of the springtails. The model proceeds in daily time steps comprising the following processes: updating the season, foraging including avoidance of contaminated cells and density dependence (hourly time steps), re-growth of food, ageing and growth, reproduction, hatching, density dependence on fecundity and survival, and mortality.
Values of almost all parameters are drawn from uniform or normal probability distributions, in order to reflect heterogeneity among individuals. Stochasticity is also used for initializing springtails' starting positions, as well as causing individual behaviours (movement, reproduction, hatching, mortality) to occur with specified frequencies, which depend on the values of said parameters.
Food resources are also randomly assigned at the beginning of a model run to grid cells which are initialised to be food sources, with different maximal food levels.
A full description of the model following the ODD format (Grimm et al., 2006; Grimm et al., 2010) is provided in the Supplementary Material.
Submodel generation of fractal patterns of habitat destruction
This submodel is based on a NetLogo implementation of the midpoint displacement algorithm (Saupe, 1988) included in the individual-based model "TraitScape" developed by Jackson and Fahrig (2012), which has been modified to fit our purpose.
The midpoint displacement is a well-known algorithm that produces random, realistic-looking, fractal landscapes (Saupe, 1988). The amount of spatial autocorrelation of the fractals generated by this algorithm can be adjusted by varying the value of a parameter (H). The fractal dimension (D) of the landscape is a property of H such that D = 3 - H (Jackson and Fahrig, 2012). Another characteristic that makes this algorithm especially suitable to study natural phenomena is that it allows independent control of habitat amount and configuration. The main differences among landscapes in our runs are the amount and the configuration of habitat, which are driven by the following parameters:
Habitat cover. Proportion of landscape cells that are clean habitat. The user can control habitat cover by adjusting a parameter (user-cover) to the desired amount. For instance, a user-cover value of 0.3 means that 30% of the habitat is without contamination.
Spatial autocorrelation of habitat (H). Degree of aggregation of habitat cells, i.e. the opposite of habitat fragmentation. H can assume values between 0 and 1; given the same habitat cover, low H will result in many small patches and low inter-patch distances, whereas high H will result in a few large patches with high average inter-patch distances (Fig. 1a-c).
H=0 color1.gifH=0.72 color1.gifH=1 color1.gif
Figure a-c Examples of fractal landscapes with 20% habitat cover (dark grey), and different degrees of spatial autocorrelation of the clean habitat: uncorrelated (a), moderately correlated (b) and completely correlated (c)
Submodel disturbances
Disturbance events are characterized by disturbance level, number of events in a year, days of occurrence of disturbance events, and spatial autocorrelation of disturbed area. The disturbance level is defined as the proportion of area disturbed per single disturbance event. The shape and placement of disturbed areas are either randomized, in the case of uncorrelated disturbances, or distributed in a circle, of area proportional to the disturbance level and located around the central grid cell, in the case of spatially autocorrelated disturbances. Disturbances occur with a frequency defined by the specified number of events, removing all individuals (eggs, juveniles and adults) on the disturbed grid cells, and can hit both habitat cells and unsuitable (i.e. contaminated by copper) cells.
Submodel summer drought
Implementation of the summer drought submodel is based on data from Waagner et al (2011). In this study, the authors performed a long term experiment which the water potential of soil was slowly decreased, to reproduce the natural condition that occurs when soil dries out. Relative humidity (RH) was progressively decreased during 12 days and the target level was maintained for 20 days. Exposure to RH > 98.2% had no significant effect on survival, while below 99.4% RH oviposition stopped. Among the different target RH tested in this study, 97% RH has been chosen for implementation, as the range 99.8 to 97% RH represents a realistic RH regime in soil during periods of natural drought (Holmstrup, 1997; Hojer et al., 2001). Furthermore, data from a study by Holmstrup (1997) showed a reduction in drought tolerance caused by copper when the desiccation stress is higher than 97.8% RH.
Therefore, based on these observations, the implemented drought effects reflect the decline in relative humidity shown in Fig. 2, and the following assumptions have been made:
Four days after the beginning of the drought period, correspondent to a decline of RH below 99.4%, all eggs die and both juvenile and adult survival begin to be affected by drought.
As reported by Waagner et al. (2011), mean survival at the end of the exposure period to 97% RH is 32%. Survival is therefore implemented as the probability to survive each day until the end of the drought period (i.e. 25 days, from day 5 to day 29):
To account for the variability in drought tolerance, individual survival follows a normal distribution with the same mean and standard deviation as recorded in Waagner et al. (2011).
Starting from the seventh day of drought, equivalent to a RH value of 97%, drought tolerance is reduced by 30% if the individual has cumulatively been exposed to copper for at least a week, according to Holmstrup (1997). This rule does not apply to the first six days of drought, as no reduction of drought tolerance was observed for RH values above 97.8% (Holmstrup, 1997).
