As a result of the Gorkha earthquake in 2015, about 9000 people lost their lives and many more were injured. Most of these losses were caused by earthquake-induced landslides. Sustainable planning and decision-making are required to reduce the losses caused by earthquakes and related hazards. The use of remote sensing and geographic information systems (GIS) for landslide susceptibility mapping can help planning authorities to prepare for and mitigate the consequences of future hazards. In this study, we developed landslide susceptibility maps using GIS-based statistical models at the regional level in central Nepal. Our study area included the districts affected by landslides after the Gorkha earthquake and its aftershocks. We used the 23,439 landslide locations obtained from high-resolution satellite imagery to evaluate the differences in landslide susceptibility using analytical hierarchy process (AHP), frequency ratio (FR) and hybrid spatial multi-criteria evaluation (SMCE) models. The nine landslide conditioning factors of lithology, land cover, precipitation, slope, aspect, elevation, distance to roads, distance to drainage and distance to faults were used as the input data for the applied landslide susceptibility mapping (LSM) models. The spatial correlation of landslides and these factors were identified using GIS-based statistical models. We divided the inventory into data used for training the statistical models (70%) and data used for validation (30%). Receiver operating characteristics (ROC) and the relative landslide density index (R-index) were used to validate the results. The area under the curve (AUC) values obtained from the ROC approach for AHP, FR and hybrid SMCE were 0.902, 0.905 and 0.91, respectively. The index of relative landslide density, R-index, values in sample datasets of AHP, FR and hybrid SMCE maps were 53%, 58% and 59% for the very high hazard classes. The final susceptibility results will be beneficial for regional planning and sustainable hazard mitigation.