Forest soil organic carbon (SOC) accounts for a large portion of global soil carbon stocks. Accurately mapping forest SOC stocks is a necessity for quantifying forest carbon cycling and forest soil sustainable management. In this study, we used a boosted regression trees (BRT) model to predict the spatial distribution of SOC stocks during two time periods (1990 and 2015) and calculated their spatiotemporal changes during 25 years in Liaoning Province, China. A total of 367 (1990) and 539 (2015) sampling sites and 9 environmental variables (climate, topography, remote sensing) were used in the BRT model. The ten-fold cross-validation technique was used to evaluate the prediction performance and uncertainty of the BRT model in two periods. It was found that the BRT model could account for 65% and 59% of SOC stocks, respectively for the two periods. MAP and NDVI were the main environmental variables controlling the spatial variability of SOC stocks. Over the 25-year period, the average SOC stocks increased from 5.66 to 6.61 kg m −2 . In the whole study area, the SOC stocks were the highest in the northeast, followed by the southwest, and the lowest in the middle of the spatial distribution pattern in the two periods. Our accurate mapping of SOC stocks, their spatial distribution characteristics, influencing factors, and main controlling factors in forest areas will assist soil management and help assess environmental changes in the region. soil properties are affected by many environmental factors, it is difficult to predict soil properties accurately and efficiently at the regional scale [8]. As an efficient and low-cost method, digital soil mapping (DSM) technology is able to accurately describe the spatial variability of soil attributes in a region by using a limited amount of sampling data and environmental variables [9]. In order to accurately predict the SOC stocks in a region and analyze its key environmental factors, scholars have carried out a lot of researches [10][11][12]. In fact, the direct relationship between soil properties and environmental factors is non-linear and complex [13,14]. Therefore, machine learning algorithms which can effectively avoid this kind of problem are widely used in DSM mapping [11,15]. For example, on Barro Colorado Island, Panama, Grimm et al. [15] used a random forest model to simulate SOC stocks in 0-10 cm, 10-20 cm, 20-30 cm, and 30-50 cm layers, respectively. Giasson et al. [16] applied a multiple logistic regression model to predict the occurrence of soil types in southern Brazil. Dorji et al. [17] selected regression kriging and equal-area spline profile function to predict the SOC stocks at depths 0-5, 5-15, 15-30, 30-60, and 60-100 cm in montane ecosystems, Eastern Himalayas. In Denmark, Adhikari et al. [18] used regression kriging and 12 environmental variables to simulate the spatial distribution of SOC stocks at five soil layers. Were et al. [19] compared the performance of support vector regression, artificial neural network, and random forest models in predicti...