Nematodes are indicators of soil quality and soil health. Knowledge of the relationships between nematode-based soil quality indices and environmental properties is beneficial 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 significant portion of the spatial variation of the indices and that regression kriging could more accurately model their spatial distribution than ordinary kriging. A stratified 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 fitted 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 significant 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.