The multilayered soil slopes are very common in geotechnical practice. Not only there exists inherent spatial variation for soil properties, but also there shows obvious multilayered distribution characteristics within the natural soil profile. However, few attempts have been made to study the reliability of multilayered soil slope stability in spatially variable soils. This paper proposes a Monte-Carlo simulation (MCS) based slope system reliability analysis approach using multiple response surfaces. The flow chart for the proposed approach is presented. The reliability of a multilayered soil slope system considering spatial variability of soil properties is systematically investigated using the proposed approach. The results indicate that the proposed approach properly evaluates the slope system reliability at small probability levels (i.e., pf,s < 0.001) considering spatial variability of soil properties and achieves high computational efficiency in parametric sensitivity analysis.