Abstract The Maharloo watershed has witnessed many gullies in the recent due to the specific topo-climatic conditions and man-made activities in that area. The present study is set out to address this issue by producing gully erosion prediction maps via three machine learning algorithms including RF, SVM and BRT in Maharloo watershed, Fars province, Iran. Also, this research attempted to consider the importance of effective factors in the occurrence of gully erosion using Boruta algorithm. To this end, gully erosion locations were identified by extensive field surveys as well as the use of already prepared gully raster map of Maharloo watershed. Then, sixteen causative factors of gully erosion such as elevation, slope degree, slope aspect, plan curvature, TWI, distance from rivers, distance from roads, drainage density, lithology, annual mean rainfall, NDVI, land use and some soil characteristics (pH, clay percent, electrical conductivity-EC, and silt percent) were identified and their maps were produced and classified in the GIS. In this study, the relationships among each agent and gully erosion were defined employing the evidential belief function (EBF) algorithm and the weight of each factor's classes was determined. On the other hand, the results of the collinearity test among the factors showed that sand percentage agent had a VIF > 5; therefore, this covariate was removed from the model. Also, the results of the importance of effective factors using Boruta algorithm indicated that three factors including land use, distance from river, and clay percent had the most noticeable importance in the occurrence of gully erosion in the study area. Finally, the gully erosion susceptibility maps were produced using the RF, BRT, and SVM models in the R statistical software. The results of machine learning techniques were evaluated employing 30% of unused locations in the modeling process as well as the receiver operating characteristic (ROC) curve. Also, in the current research, try to assess the fitting performance of models and their robustness using sensitivity rate, specificity rate, Cohen's Kappa, and 4-fold cross-validation measures. Results showed that the final gully erosion susceptibility maps had an excellent accuracy with AUC values of validation data sets by scenarios 7 and 9 independent factors on gully erosion, respectively (SVM = 0.957, 0.975, RF = 0.991, 0.986, BRT = 0.913, 0.913). The fitting performance measures and robustness technique (4-fold cross-validation) also confirmed the achieved validation results. In order to control and prevent this type of erosion in the Maharloo watershed, there should be protective actions and watershed management measures in place at the primary stages, especially at the beginning of the gully erosion, to control the development of the gully erosion.