Numerical simulation for large deformation problems of saturated soil poses an ongoing challenge since it requires the capabilities of modeling large deformations, interactions between solid skeleton and pore water, as well as fluidization behaviors. A hydro-mechanical coupled B-spline material point method is developed to simulate large deformations of saturated soils. Based on Biot's mixture theory, a single set of material particles is adopted to represent saturated soils where Darcy's law is employed to govern the motion of pore water and the momentum conservation of the mixture is used to govern the motion of solid skeleton. The convergence and effectiveness of the proposed method are examined by one-dimensional consolidation tests and quasi-static footing penetration into two-dimensional elastic saturated soils. The capability of simulating practical engineering for the proposed method is demonstrated by reproducing the dynamic process of three-dimensional saturated soil slopes with a strain-softening Drucker-Prager constitutive model. Results show that the coupled BSMPM well reveals the overall dynamic process of saturated soil slope involving evolution of shear bands, formation of graben and horst, and fluidization of progressive failures with improved computational accuracy and enhanced convergence, which demonstrates the proposed method is a promising approach for large deformation problems of saturated soils.