Geostatistical modelling and mapping of nematode-based soil ecological quality i

Israel O.IKOYI*Gerard B.M.HEUVELINK and Ron G.M.DE GOEDE

1 Wageningen University&Research，Soil Biology Group，Wageningen 6700 AA(The Netherlands)

2 Current address:University of Limerick，School of Natural Sciences，Department of Biological Sciences，Limerick V94T9PX(Ireland)

3 Wageningen University&Research，Soil Geography and Landscape Group，Wageningen 6700 AA(The Netherlands)

ABSTRACT Nematodes are indicators of soil quality and soil health.Knowledge of the relationships between nematode-based soil quality indices and environmental properties is benef icial for assessing environmental threats on soil biota.This study evaluated the spatial distribution of nematode-based soil quality indices in a 23-ha heavy metal-polluted nature reserve using geostatistical methods.We expected that a selection of abiotic soil properties(pH and moisture，clay，organic matter，cadmium(Cd)，and zinc(Zn)contents)could explain a signif icant portion of the spatial variation of the indices and that regression kriging could more accurately model their spatial distribution than ordinary kriging.A stratif ied simple random sampling scheme was used to select 80 locations where soil samples were taken to extract nematodes and derive the indices.The area had a distinct gradient in soil properties with Cd and Zn content ranging from 0.07 to 68.9 and 5.3 to 1 329 mg kg-1，respectively.Linear regression models were f itted to describe the relationships between the indices and soil properties.By also modelling the spatial correlation structure of regression residuals using spherical semivariograms，regression kriging was used to produce maps of the indices.The regression models explained between 21%and 44%of the total original variance in the indices.Soil pH was a signif icant explanatory variable in almost all cases，while heavy metal conent had a remarkably low effect.In some cases，the regression residuals had spatial structure.Independent validation indicated that in all cases，regression kriging performed slightly better because of having lower values of the root mean square prediction error and a mean prediction error closer to zero than ordinary kriging.This study showed the importance of soil properties in explaining the spatial distribution of biological soil quality indices in ecological risk assessment.

Key Words: ecological risk assessment，heavy metals，model validation，regression kriging，semivariance analysis，soil property，spatial structure

INTRODUCTION

In terrestrial ecosystems，soil biodiversity consists mainly of microorganisms and small invertebrates(Lavelle and Spain，2001).Studies have shown that the distribution of these organisms is not random，but displays spatial patterns over scales ranging from square millimeters to hectares(Görreset al.，1998;Ettemaet al.，2000;Stoyanet al.，2000;Berg，2012;Frey，2015;Songet al.，2017).Competition，habitat conf iguration，and historical effects have been identif ied as possible factors inf luencing the spatial distribution of soil organisms and spatial structure of soil biota communities(Martinyet al.，2006;Horner-Devineet al.，2007;Vos and Velicer，2008).Analysis of soil biota community structure can provide valuable information on ecosystem health(Kibblewhiteet al.，2007;Morales-Márquezet al.，2018).However，the spatial complexity of soil biota communities complicates adequate determination of bioindicators and subsequent ecosystem health assessment.Ecosystem health assessment studies often usead hocsampling designs，where decisions on the number of soil cores per sampled surface are arbitrary(Quistet al.，2017)and widely spaced soil samples are bulked into composite samples.A consequence of the latter is that biota assemblages，rather than communities of specimens that can potentially interact，are studied，which questions the use of such data in food web relation studies.Hence，scientif ically well-considered sampling strategies are needed to improve the quality of the data that are used to provide spatially explicit soil quality assessment indicators.

Nematodes have been reported to be the most abundant metazoans on the earth’s surface with about 1×108individuals distributed in a square meter of soil(Lambshead，2004;Decraemer and Hunt，2006).They are assumed to occur in almost all environments，in nearly every type of soil，under all climatic conditions and in all ecological niches，varying from undisturbed to disturbed ones(Bongers and Ferris，1999).Different nematode taxa utilize specif ic food sources，and shifts in these feeding groups often ref lect changes in the soil food web(Ferriset al.，2001;Yeateset al.，2009).In addition，nematode species differ in their sensitivities to different forms of stress and disturbance，including toxic effects due to heavy metal contamination(Korthalset al.，1996;Georgievaet al.，2002;Tenuta and Ferris，2004;Sánchez-Moreno and Navas，2007).The unique combination of feeding habit，life span，and sensitivity to environmental disturbance has allowed the classif ication of nematode taxa with similar response characteristics into functional guilds(Bongers and Bongers，1998;Ferriset al.，2001;Ferris and Bongers，2006).This was used to develop nematode-based indices which are deployed as bioindicators of environmental health，soil quality，and ecosystem resilience(Bongers，1990;Ferriset al.，2001;Neher，2001;Höss and Traunspurger，2003;Yeateset al.，2003;Mulderet al.，2005;Schratzbergeret al.，2006;Heiningeret al.，2007).

The spatial distribution of nematodes can be related to spatial heterogeneity in soil properties(Robertson and Freckman，1995;Ettemaet al.，1998;Lianget al.，2005;Monroyet al.，2012;Park，2012;Quistet al.，2017;Songet al.，2017)and is sometimes inf luenced by differences in land use intensity and pollution.Improvement of our understanding of spatially dependent relationships between nematode-based soil quality indices and environmental variables can contribute to the practical application of such indices in studies on(precision)agriculture，landscape ecology，biodiversity，etc.Within an agricultural and ecological setting，it is important to know the consequences of the tradeoffs between plant parasitic nematodes and nontarget organisms.

Geostatistics provides powerful tools with its basic functions，such as the variogram，used for spatial analysis and interpolation.It quantif ies spatial correlation，which can include the inf luence of explanatory variables in a trend，interpolates spatial data using kriging(Cressie，1990;Oliver and Webster，2014;Rosemaryet al.，2017)，and quantif ies the interpolation error with the kriging variance(Oliver and Webster，2014).Most soil biological studies that utilize geostatistical tools have focused on specif ic nematode genera or feeding groups(Wallace and Hawkins，1994;Lianget al.，2005;Ortizet al.，2010).Merckxet al.(2010)used geostatistics to study free-living marine nematodes in the Southern Bight of the North Sea and produced biodiversity maps;however，their study focused on diversity indices and did not consider the inf luence of disturbances，such as heavy metal contamination.Extending such studies on nematodebased soil quality indices in an area under disturbance is valuable as it can reveal the inf luence of the disturbance on the nematode community composition and improve our understanding of the mechanisms that cause spatial variation within the community.With this background，there is a need to study the spatial distribution of nematode-based indices in relation to environmental properties(including human heavy metal contamination).With a strategic spatial sampling design to collect a limited number of samples，the impact of environmental pollution can be assessed more efficiently.This study was designed to achieve the following objectives:i)describe the spatial structure of nematodebased indices in a natural area with different levels of soil pollution;ii)establish relationships between these indicators and soil properties including heavy metal content;iii)use these relationships to construct spatially explicit maps of the indices in a regression kriging approach;and iv)validate and interpret the resulting models and maps.Geostatistical methods were applied to achieve these objectives using a unique dataset from an area with distinct gradients in abiotic soil conditions and pollution levels.

MATERIALS AND METHODS

Study area

This study was conducted on a 23-ha area(5°27′E，51°18′N)in the Malpiebeemden Nature Reserve，located in the province of North-Brabant in the south of The Netherlands(in between the municipality of Valkenswaard and the Dutch-Belgian border)(Fig.1a).The study area comprised grassland，reed land，and some small forest patches west of the Dommel River.Due to upstream mining activities，the Dommel River was polluted with heavy metals such as cadmium(Cd)and zinc(Zn)until recently(Bleeker and van Gestel，2007).Regular f looding events in the past resulted in increased heavy metal concentrations in large parts of the area.The highest topsoil Cd concentrations(up to 100 mg kg-1)in The Netherlands are found in the study area.The distribution of the contamination in the f loodplains is very heterogeneous，resulting in large gradients of heavy metal concentrations(Bleeker and van Gestel，2007).This area was selected because of the large heterogeneity in soil pollution levels and also because it exhibits large differences in basic soil properties，such as soil texture，organic matter，pH，and moisture.The vegetation of the nonforested part is dominated byHolcuslanatus，Juncuseffusus，andPhragmitesaustralis.There are four different soil types(Albic Arenosol，Gleyic Podzol，Humic Podzol，and Peat soils)in the study area.

Soil sampling

Soil samples for nematode analyses were collected on March 20，2013 using a stratif ied simple random sampling design.This sampling design was used because it guarantees a fairly uniform spread of sampling locations over the area and provides data suitable for model validation(de Gruijteretal.，2006).It also yields short-distance comparisons necessary for semivariogram f itting.The study area was divided into eight strata of equal size using the spcosa package in R(de Gruijteret al.，2006;Walvoortet al.，2010).Prior to stratif ication，the forested patches(Fig.1b)were masked in order to have soil samples from only uniform vegetation(i.e.，grasslandsensu lato).Ten soil sampling points were chosen randomly within each stratum，making a total of 80 soil sampling points(Fig.1b).Moreover，before visiting the f ield，f ive extra random soil sampling points per stratum were selected as a precaution in case any of the preselected main sampling points within the stratum turned out to be in an unsuitable location(e.g.，inside a forest patch or pond).These f ive additional soil sampling points had a ranking order，and the rank was taken into consideration when any of the main sampling points had to be replaced，which happened in three cases.In the f ield，the positions of the soil sampling points were located with a geographical positioning system(GPS).At each of the sampling locations，one soil sample was taken from a soil depth of 0—20 cm using a soil corer with a diameter of 4 cm.At each date of sampling，the air temperature was close to 0°C.The samples were stored at 4°C until nematode extraction.

Soil samples for physico-chemical analyses(pH，clay，organic matter，moisture，and reactive and soluble Cd and Zn)were collected on April 16 and 18，2008.Using a spatial grid with a internode distance of 50 m，a total of 77 soil samples at 0—20 cm soil depth were collected.In addition，another 23 soil samples were collected at internode distances of 25 and 12.5 m.The samples were dried at 40°C and stored for further analysis.Soil pH was measured in a 0.01 mol L-1CaCl2solution(1:10，weight/volume).Soil organic matter content was measured using loss-on-ignition at 550°C.Soil moisture content was determined using oven drying at 105°C overnight.Soil clay content was measured by BLGG AgroXpertus(Wageningen，The Netherlands).Soil reactive Cd and Zn contents were measured with a f lame atomic absorption spectrometer(AAS)(AAnalist 300，Perklin Elmer，Norwalk，USA)，after extraction in 0.43 mol L-1nitric acid solution(1:10，weight/volume).Soil soluble Cd and Zn were extracted in 0.01 mol L-1CaCl2solution(1:10，weight/volume)and measured with the f lame AAS or in case of low contents with an inductively coupled plasma mass spectrometer(ICP-MS，Element 2，Thermo Scientif ic，Waltham，USA).

Nematode extraction and identif ication

Nematodes were extracted from 100 g fresh soil samples using the Oostenbrink elutriator(Oostenbrink，1960)within three weeks after sampling.Nematodes from 10%of the total extracted volume were counted using a high magnif ication microscope at 100—400 times magnif ication.The total number of nematodes was expressed as individuals kg-1dry soil，taking the soil moisture content into account.They were heat killed and f ixed in 4%formaldehyde and identif ied to the genus or family level based on morphological features using identif ication keys(Siddiqi，1986;Bongers，1988;Jairajpuri and Ahmad，1992).This was the highest achievable taxonomic resolution and is typically the f inest resolution used in studies that report composition of soil nematode communities and is sufficient for calculation of nematode-based soil quality indices.

Nematode-based soil quality indicators

The nematodes were allocated to feeding groups according to Yeateset al.(1993)and colonizer-persister(c-p)groups using the method described in Bongers and Bongers(1998).The total nematode abundance(individuals kg-1soil)，abundance of feeding groups(individuals kg-1soil)，genera richness，Shannon index，and the maturity，enrichment，channel，and structure indices were calculated for each of the 80 samples collected on April 16 and 18，2008.

The maturity index was calculated for the non-plantfeeding taxa as:

where MI is the maturity index，viis the c-p value assigned to familyi，f iis the frequency of familyi，andnis the total number of individuals of the non-plant-feeding taxa in the sample(Bongers，1990).

Structure，enrichment，and channel indices were calculated following Ferriset al.(2001)and Ferris and Matute(2003).The basal(b)，enrichment(e)，and structure(s)components of the nematode assemblages were calculated as:

where BA1，BA2，and BAmrepresent the abundances of bacterial-feeding nematodes belonging to c-p 1，2，and m，respectively，FU2and FUmand CA1and CAmrepresent the abundances of fungal-feeding and carnivorous nematodes belonging to c-p 2 and m，respectively，OAmrepresents the abundance of omnivorous nematodes belonging to c-p m，andW1，W2，andWmare the weights assigned to nematodes belonging to c-p 1，2，and m，respectively(Ferriset al.，2001).Subsequently，the enrichment and structure indices were calculated as:

where EI is the enrichment index and SI is the structure index.The channel index is the proportion of fungal-feeding nematodes in c-p 2(FU2)within the opportunistic decomposer guilds，including fungal-feeding nematodes in c-p 2 and bacterial-feeding nematodes in c-p 1(BA1)(Ferriset al.，2001).Taking into consideration the specif ic weights of each guild，it was calculated as:

where CI is the enrichment index.The genera richness was calculated as the number of genera in each sample，and the Shannon index(H)was calculated as:

wherePiis the proportion of generaiin the sample.

Semivariogram and kriging

Before conducting the geostatistical analysis，an exploratory data analysis was carried out by computing standard summary statistics and graphics.The spatial dependence of abiotic soil properties and indices and associated regression residuals were quantif ied using the sample semivariogram，ˆγ(h):

whereN(h)is the number of paired observations separated by distanceh，z(si)is the value of the variable of interest at locationsi，andz(si+h)is the value of the variable of interest at a location at distancehfromsi(Webster and Oliver，2007).Next，semivariogram models were f itted to the sample semivariograms using weighted least squares.A semivariogram model typically has three important parameters，sill，nugget，and range.The sill is the maximal semivariogram value and a measure of the total spatial variation(and similar to the variance of the data).The nugget is the combined effect of random measurement errors and short-distance spatial variation(Goovaerts，1999;Guanet al.，2017;Yaoet al.，2019).The range is the distance where the semivariogram model reaches the sill，or the distance up to which there is spatial correlation.The relative structure of the semivariogram models，i.e.，the nugget-to-sill ratio，was calculated for each of the nematode-based indices(Robertson and Freckman，1995).Nematode data that showed no spatial structure(i.e.，pure nugget effect，based on visual observation of the sample semivariogram and f itted parameters)were excluded from further analyses(modelling and mapping).

Geostatistical tools in R(R Core Team，2018)were used to calculate and f it the semivariograms of the residuals from the regression model to test for spatial autocorrelation.The relative structure(nugget-to-sill ratio)of the semivariograms(i.e.，the proportion of the sample variance that is spatially autocorrelated)was calculated for each of the residuals that showed spatial autocorrelation(Robertson and Freckman，1995).

Ordinary kriging was used to interpolate the abiotic soil parameters measured in 2008 at 100 sampling points to the 80 sampling points measured in 2013.Soil clay，organic matter，and heavy metal contents had hardly changed during the 5-year intervening period.This may，however，not be true for soil pH and moisture content.However，since abiotic soil properties were used as explanatory variables in the regression models，it is only their relative spatial differences that are important，not their absolute values.These spatial differences will stabilize over time.In addition，there is no theoretical objection against using explanatory variables from f ive years earlier to predict a dependent variable.Perhaps the strength of the relationship between the dependent and independent variables is weakened by the time difference—although there are cases where the effect of environmental change on organisms is better modelled with a time-lagged regression approach(e.g.，Ludovisiet al.，2014);otherwise，the regression is as valid as when all variables were measured the same time.

Modelling relationships between the nematode-based indices and abiotic soil properties

Multiple linear regression models were used to establish the relationships between the indices，including maturity index，structure index，Shannon index，and genera richness，as dependent variables，and soil properties，including pH and organic matter，moisture，clay，Cd，and Zn contents，as explanatory variables.As both reactive and available Cd and Zn contents showed high positive correlation(r2＞0.9)，only the reactive contents were used in model building.The explanatory variables were included using forward stepwise selection based on the Akaike information criterion(AIC).The general form of the multiple linear regression models was:

whereYis the dependent variable，β0is the intercept，β1—βpare the coefficients of the independent variablesX1—Xp，andεis the stochastic residual，which is assumed to be normally distributed with zero mean and constant variance.

Negative binomial regression models were used for modelling the relationships between countable dependent variables(total abundance and abundance of feeding groups)and soil properties(Zeileiset al.，2008)，using the glm.nb()function in the MASS package in R(Venables and Ripley，2002).In the negative binomial regressions，the dependent variableYhad a negative binomial distribution with its mean(μ)given by:

Geostatistical interpolation of the nematode-based indices was done using a regression kriging approach.The generic model for regression kriging was also given by Eq.10，but unlike in multiple linear regression，in this case，the residualεwas checked for possible spatial dependence.If the residuals were spatially dependent，the semivariogram was estimated，simple kriging of residuals(assuming zero mean)was used to interpolate the residuals，and the interpolated residuals were added to the regression prediction(Henglet al.，2004).This was done to improve the prediction of the dependent variable.

Correlation analysis

Pearson correlation coefficient values were calculated in R(RCore Team，2018)to examine the relationships between soil properties and the nematode-based soil quality indices.Correlations atP＜0.05 were considered as signif icant.

Validation of the per formance of prediction maps

The performance of the prediction maps was validated using leave-one-out cross-validation(Webster and Oliver，2007).Prediction was evaluated by comparing the predictions with actual observations at validation points to assess the systematic error，calculated as mean prediction error(MPE)，and the prediction accuracy，calculated as the root mean square prediction error(RMSPE)(Henglet al.，2004).We also computed the normalized RMSPE，def ined as the square root of the average ratio of the squared prediction error to the kriging variance.Ideally，it should be close to one，indicating that the kriging standard deviation is a proper measure of the interpolation error.In an ideal scenario，exact interpolation is indicated by a zero kriging standard deviation(Oliver and Webster，2014)，which is，however，unlikely in real-world applications(Hatvaniet al.，2017).The validation metrics were compared with those obtained when using ordinary kriging of the indices，in order to analyse the added value of including information derived from abiotic soil properties in mapping.

RESULTS

Soil physicochemical properties and nematode assemblages

The soils tested were acidic，with pH values ranging from 3.7 to 4.9(Table I).The clay，moisture，and organic matter contents were highly variable.The mean values of reactive Cd(10.3 mg kg-1)and Zn(295.1 mg kg-1)contents were above the Dutch target values for these metals(0.8 and 140 mg kg-1，respectively).The highest values of most of the soil properties were recorded in the eastern part of the study area，at locations closer to the river(Fig.2).

Fig.2 Interpolated maps of abiotic soil properties:pH(a)，clay content(g kg-1，b)，organic matter content(g kg-1，c)，moisture content(g kg-1，d)，reactive Cd content(mg kg-1，e)，and reactive Zn content(×102 mg kg-1，f).

A total of 69 nematode genera belonging to 33 families were identif ied.The total number of nematodes varied from 750 to about 100 000 individuals kg-1dry soil(Table II).The most prevalent genera werePlectus，Aphelenchoides，andEucephalobus(occurring in＞90%of the samples).The maturity index values ranged from 1.3 to 2.8，the Shannon index from 1.1 to 3.0，and the genera richness from 8 to 32.The values of structure，enrichment，and channel indices，ranged from 0 to 88，4 to 93，and 0 to 100，respectively.Bacterial feeders were the most abundant trophic group(about 50%of the total nematode abundance)，whereas omnivores and predators were the least abundant(about 3%).

TABLE I Descriptive statistics of physical and chemical properties of the soils sampled on March 16 and April 18，2008 in the study area(n=100)

Semivariograms of nematode-based indices and abiotic soil data

All nematode-based indices except the enrichment and channel indices were spatially structured(Table III).The semivariogram ranges of the nematode-based indices that showed spatial dependence varied at the range values of 80—161 m.The semivariograms of maturity and structure indices and their residuals are shown in Fig.3;the results for other spatially structured indices were similar and therefore are not shown.All abiotic soil properties were spatially structured at the range values of 57—249 m(data not shown).

Fig.3 Semivariograms of the maturity index(a)and structure index(b)of nematodes and of the residuals from the regression models for the maturity index(c)and structure index(d)for the soils sampled on March 20，2013 in the study area.

Regression models

No statistically signif icant regression models(P＞0.05)were obtained for the enrichment and channel indices.The regression models for the remaining indices were statistically signif icant(P＜0.001).The multiple linear regression models f itted to the maturity index，structure index，Shannon index，and genera richness explained 21%—27%of the variation in these indices(Table IV).Soil pH was a statistically signif icant explanatory variable in all four models，and the index values decreased with increasing soil pH.Negative binomial regression models were used to model the variation in the remaining indices and reduced the null deviance by 26%—44%.In this case，soil pH was included in f ive out of six models，again showing a negative correlation with the indices.Overall，soil pH was selected as a signif icant explanatory variable in nine out of ten models and always had a negative coefficient value.Soil heavy metal contents were selected in only four of the ten models and had variable correlation coefficient values.

The residuals of eight out of ten regression models showed spatial structure.Between 33%and 80%of the residual variance was spatially dependent at the semivariogram range values of 28—112 m(Table V).

TABLE II Descriptive statistics of nematode-based indices for the soils sampled on March 20，2013 in the study area(n=80)

TABLE III Semivariogram characteristics of the nematode-based indices for the soils sampled on March 20，2013 in the study area

Regression kriging and associated standard deviation maps，model validation，and comparison analysis

The regression kriging models were used to create maps for all nematode-based soil quality indices.The results for the maturity index and structure index are shown in Fig.4;the results were similar for other indices and therefore are not shown.The lowest maturity index predictions(＜2.0)were obtained at locations close to the river，while maturity index predictions higher than 2.4 were obtained at locations farther from the river and near the northern border of the area(Fig.4b，e).A similar pattern was observed for the structure index，where predictions higher than 60%were obtained at locations farther away from the river(Fig.4h，k).Using the spatial components in the regression model residuals(Fig.4d，j)for a further improvement of the maps(Fig.4e，k)resulted in a more detailed map of the spatial distribution of both indices.The kriging standard deviation maps(Fig.4f，l)showed that the uncertainty associated with the predictions ranged from10%to 11%for the maturity index and from 31%to 46%for the structure index.The highest accuracy levels were achieved in areas with a relatively high observation density.

Fig.4 Maps of maturity index(a—f)and structure index(g—l)of nematodes for the soils sampled on March 20，2013 in the study area:values at observation points(a and g)，ordinary kriging predictions(b and h)，regression-only predictions(c and i)，interpolated regression residuals(d and j)，regression kriging predictions(e and k)，and regression kriging standard deviations(f and l).

TABLE IV Summary of regression analyses for the nematode-based indices with abiotic properties for the soils sampled on March 20，2013 in the study area

TABLE V Semivariogram characteristics of the regression model residuals showing spatial autocorrelation for the nematode-based indices for the soils sampled on March 20，2013 in the study area

In general，the values of RMSPE and normalized RMSPE were relatively low for regression kriging compared to ordinary kriging(Table VI).For instance，the value of RMSPE for the maturity index obtained with regression kriging was 0.26 compared to 0.29 obtained with ordinary kriging.

DISCUSSION

The present study focused on the spatial distribution of nematode-based soil quality indices and their dependence on and relationship with abiotic soil properties including heavy metal contents.The nematode-based indices evaluated in this study exhibited spatial patterns at the semivariogram range values of 80—161 m(Table III).These patterns of nematodebased soil quality index distribution could be attributed to patterns of their life history characteristics and responses to soil properties(Table IV).This is in agreement with the f indings of Monroyet al.(2012)，who found that spatial clustering is a characteristic feature of both bacteria and nematode communities.Other studies also showed that at various sampling scales，nematode populations are spatially patterned(Ettemaet al.，1998;Lianget al.，2005;Park，2012;Quistet al.，2017).Even in agricultural soils，which are considered relatively homogenous，spatial patterns are detected in the distribution of nematode groups(Robertson and Freckman，1995).The semivariograms of the enrichment and channel indices showed a pure nugget effect，indicating that the distribution of these indices was random and not spatially structured at the scales studied.Both indices also did not show any signif icant relationships with the abiotic soil properties.The enrichment and channel indices were calculated based on the proportions of c-p 1 bacterial feeders and c-p 2 fungal feeders that would increase rapidly upon increase in microbial activity as a result of organic input(Ferriset al.，2001).The study area was a nature reserve which was not fertilized，but young cows were used for extensive grazing.The grazing management could have resulted in high variability in the enrichment and channel indices due to unpredictable fecal inputs rather than siterelated soil characteristics.

TABLE VI Cross-validation results of the regression kriging(RK)and ordinary kriging(OK)maps of the nematode-based indices for the soils sampled on March 20，2013 in the study area

In the present study，about 21%—44%of the variation in the indices was explained by the soil properties.Although different explanatory variables tended to explain different indices(Table V)，soil pH was almost always a signif icant explanatory variable.In most of the models，soil moisture content was also a signif icant explanatory variable.Both soil pH and soil moisture content had negative relationships with the indices.The heavy metals were selected as explanatory variables in only 4 out of 10 models.This suggested that heavy metal contents had remarkably small effects on the nematode-based indices.However，some care should be taken when interpreting these results，since cross-correlations between soil properties could have led to the selection of other cross-correlated factors.For instance，the inf luences of heavy metals may be represented by soil pH and moisture content if these soil properties are strongly correlated with heavy metal contents(Sochováet al.，2006).The values of Pearson correlation coefficient of abiotic soil properties were all in the range of 0.54—0.99，lower than those of the heavy metal contents.The values of Pearson correlation coefficient between soil pH and heavy metal contents ranged from 0.54 to 0.60;thus，we cannot exclude a confounding effect.

Pearson correlation analysis indeed showed that the total nematode abundance，maturity index，structure index，Shannon index，genera richness，and abundances of plant feeders，omnivores，and predators were signif icantly negatively correlated with all abiotic soil properties.Negative correlations were also expected for the heavy metals as the toxic effects of these metals on total nematode abundance are often reported in the literature(Korthalset al.，1996;Yeateset al.，2003;Nagyet al.，2004;Zhanget al.，2007;Shaoet al.，2008).Moreover，signif icant negative correlations of maturity index，trophic diversity，and proportions of bacterial and fungal feeders with clay，organic matter，Cd，and Zn contents have been reported previously(van Vliet and de Goede，2008).Correlations of these nematode parameters with soil pH were，however，not reported in their study.The negative correlation with soil pHis surprising as it is contrary to the general expectation that an increase in soil pH has a positive impact on nematode parameters especially because the bioavailable fractions of heavy metals increase at low soil pH(Jamaliet al.，2009).This unexpected result may perhaps be partly explained by the small range of soil pH in the study area(pH of 3.7—4.9).The negative relationship could also be due to the correlation between soil pH and other soil properties.Our study added to the mixed results on the effect of soil pH on nematode abundance.Shukurovet al.(2005)，for instance，also reported a negative correlation between total nematode abundance and soil pH，while de Goede and Bongers(1994)reported a signif icant positive correlation between soil pH and the proportion of taxa that were sensitive to environmental disturbance.

Although the regression models were all signif icant(P＜0.001)，indicating the importance of the selected explanatory variables for prediction of the nematode-based soil quality indices(dependent variable)，a large portion of the variation in the indices still remained unexplained.Examination of the semivariograms of the regression residuals showed that the residuals showed spatial dependence in some cases.The spatial autocorrelation of the residuals indicated that these indices were affected by other spatially distributed factors not measured in this study.These other factors could include other biotic properties，such as plant species identity，and intrinsic population processes，such as competition for resources and reproduction(Ettema and Wardle，2002)，but only if these factors are spatially dependent.For instance，plant species identity has been reported to inf luence soil nematode assemblage composition(De Deynet al.，2004;Viketoftet al.，2005).

Another objective of this study was to include the relationships between the nematode-based soil quality indices and soil properties in mapping the spatial distribution of the indices.This was done using regression kriging.For the maturity and structure indices，the lowest values were generally obtained at locations closest to the river(eastern part of the study area).These locations have the highest values of the soil properties(pH and organic matter，clay，moisture，Cd，and Zn contents)(Fig.2).The maturity and structure indices increased with increasing distance from the river，and their highest values were obtained at the farthest distances from the river and highest elevations in the study area.Though nematodes can survive under low oxygen conditions in saturated soils(Poinar，1983)，the population growth of many taxa may be affected by accumulation of the products of anaerobic metabolism.This might be also true for this study as most of the locations nearest to the river were in saturated conditions at soil sampling.All nematode groups were the least abundant in these locations with high soil moisture.Similar effect of high soil moisture on nematode taxa has been reported previously(Ettemaet al.，1998).Higher values of maturity and structure indices indicated that the soils were in a stable condition and not stressed(Bongers，1990;Ferriset al.，2001;Viketoftet al.，2005).Similar trends were observed for the diversity indices and abundances of the feeding groups.Regarding our expectation of increased diversity being associated with low levels of disturbance，this study suggests that the nematode diversity is affected by stress conditions like high soil moisture content(saturated condition)rather than by soil heavy metal contents.Interestingly，in 3 out of 4 models where soil moisture content was selected as a signif icant explanatory variable，it had a negative relationship with the nematode-based indices.Furthermore，the f indings of this study supported those of Bertet al.(2009)，who also did not f ind signif icant relationships between soil c-p 1 nematode groups and historical pollution.

The standard deviation maps give some level of conf idence on the accuracy of the regression kriging maps.Based on the results of the cross-validation analysis(Table VI)，slightly better relative prediction accuracy was achieved with regression kriging than with ordinary kriging in all cases.This conf irmed that the soil properties explained a signif icant portion of the variation in the nematode-based soil quality indices.

CONCLUSIONS

The results of the multiple linear regressions，semivariogram analysis，and regression kriging showed that the distribution of nematode-based soil quality indices was not random，but exhibited a spatial pattern in this study.The main driving factors were soil pH and moisture content，while soil heavy metal contents had a remarkably low effect.Despite its exploratory nature，this study offers insights into the relationships between the nematode-based indices and soil properties.Geostatistical models relating the indices to abiotic soil properties and quantifying spatial correlations between them were used in producing prediction maps of the indices.Regression kriging produced slightly more accurate maps than ordinary kriging because it could benef it from the information contained in the explanatory variables.This was shown in the lower kriging standard deviation maps and conf irmed by cross-validation statistics.

This study did not reveal any direct inf luence of heavy metals(Cd and Zn)on the nematode assemblage and spatial distribution of nematode-based soil quality indices，although the inf luences of heavy metals on these indices may have been included through other cross-correlated soil properties.However，if present，this confounding effect is not likely to be very large since the values of the Pearson correlation coefficient between soil heavy metal contents and the main explanatory variable soil pH were smaller than 0.60.Thus，if soil heavy metal pollution were the main dominant factor，it might have been included in the model selection procedure much more often.Apparently，even in this heavily polluted study area，nematode abundance was less limited by heavy metal pollution and much more by other abiotic soil properties，such as soil pH and moisture content.