Digital soil mapping has been introduced as a viable alternative to the traditional mapping methods due to being fast and cost-effective. The objective of the present study was to investigate the capability of the vegetation features and spectral indices as auxiliary variables in digital soil mapping models to predict soil properties. A region with an area of 1225 ha located in Bajgiran rangelands, Khorasan Razavi province, northeastern Iran, was chosen. A total of 137 sampling sites, each containing 3-5 plots with 10-m interval distance along a transect established based on randomized-systematic method, were investigated. In each plot, plant species names and numbers as well as vegetation cover percentage (VCP) were recorded, and finally one composite soil sample was taken from each transect at each site (137 soil samples in total). Terrain attributes were derived from a digital elevation model, different bands and spectral indices were obtained from the Landsat7 ETM+ images, and vegetation features were calculated in the plots, all of which were used as auxiliary variables to predict soil properties using artificial neural network, gene expression programming, and multivariate linear regression models. According to R RMSE and MBE values, artificial neutral network was obtained as the most accurate soil properties prediction function used in scorpan model. Vegetation features and indices were more effective than remotely sensed data and terrain attributes in predicting soil properties including calcium carbonate equivalent, clay, bulk density, total nitrogen, carbon, sand, silt, and saturated moisture capacity. It was also shown that vegetation indices including NDVI, SAVI, MSAVI, SARVI, RDVI, and DVI were more effective in estimating the majority of soil properties compared to separate bands and even some soil spectral indices.