## 1 Introduction  
Soil is a basic non-renewable resource that plays an essential role in promoting sustainable agricultural development, maintaining biosphere stability, and tackling climate change [1]. Environmental degradation and food security are prominent problems in agricultural development, which calls for the effective management and protection of soil ecosystems [2]. A watershed is an area of land that drains water into a specific water body. The watershed is considered fundamental to ecological balance, agricultural production,and climate change [3]. Soil properties, such as organic carbon (SOC), total phosphorus (STP), and available phosphorus (SAP) contents, have a major influence on soil fertility, carbon cycling, and environmental quality [4; 5]. Investigating the distribution of multiple soil properties in watersheds as hydrologic units can provide insights into the latest status of soil health, allowing for the effective evaluation of agricultural production capacity. Traditional soil mapping mainly involves manual delineation based on topographic maps and aerial images, which is labor-intensive and has a limited accuracy [6]. The rapid development of 3S technology and access to open environmental data (e.g., relief, remote sensing, and climate) have enabled a growing number of digital soil mapping studies to accurately predict the distribution patterns of soil properties on regional scales [7]. Compared to traditional methods, digital soil mapping techniques produce more objective and replicable results with fewer samples [8]. Machine learning models, exemplified by the random forest (RF) [9] and extreme gradient boosting (XGBoost) models [10], are widely used to predict soil properties. The RF model has been shown to be an effective prediction method, as it prevents over-learning and provides stable predictions [11; 12]. For example, [PERSON] et al. [13] mapped a set of soil properties (e.g., bulk density, SOC, total nitrogen, and cation exchange capacity) at six soil depths using a quartile RF model with 400 environmental covariates. [PERSON] et al. [14] leveraged the RF model with Earth observation data and environmental auxiliary information to map the SOC content and storage, pH, and cation exchange capacity in the topsoil of agricultural land in Brazil. The successful application of XGBoost in soil mapping has also been documented, given its fast computation and excellent performance [15; 16]. [PERSON] et al. [17] used 10 different methods, including XGBoost, to predict the distribution of available cadmium in soils. Among the traditional and integrated models used, the XGBoost model yielded the best prediction performance, which was superior to the linear regression model in terms of accuracy. The previous studies have demonstrated the applicability of machine learning models represented by RF and XGBoost for predicting soil properties. These models provide promising tools for the accurate mapping of the soil properties in watersheds. During digital soil mapping, the number of environmental covariates increases rapidly [18]. Therefore, variable selection is commonly used before fitting the final prediction model. This process can enable rapid model calibration, reduce model complexity, improve prediction performance, avoid multicollinearity, and facilitate map generation [19]. The recursive feature elimination (RFE) algorithm is one of the most popular variable selection methods [13; 20; 21; 22; 23; 24]. When using the RF model to assess variable importance and to predict the SOC content, [PERSON] et al. [25] selected the most accurate variable combination for each set of environmental covariates based on RFE. Furthermore, [PERSON] et al. [26] assessed the prediction accuracy of images for soil organic matter across different time periods using the RF model. They determined the performance of RF models with two different variable selection methods, and the model with RFE achieved a higher accuracy.