A first-order, Monte Carlo ensemble method has been recently introduced for solving parabolic equations with random coefficients in [Luo and Wang, SIAM T. Nurner. Anal., 56 (2018), pp. 859-876], which is a natural synthesis of the ensemble-based, Monte Carlo sampling algorithm and the ensemble-based, first-order time stepping scheme. With the introduction of an ensemble average of the diffusion function, this algorithm leads to a single discrete system with multiple right-hand sides for a group of realizations, which could be solved more efficiently than a sequence of linear systems. In this paper, we pursue in the same direction and develop a new multilevel Monte Carlo ensemble method for solving random parabolic partial differential equations. Comparing with the approach in [Luo and Wang, SIAM T. Numer. Anal., 56 (2018), pp. 859-876], this method possesses a second-order accuracy in time and further reduces the computational cost by using the multilevel Monte Carlo sampling method. Rigorous numerical analysis shows the method achieves the optimal rate of convergence. Several numerical experiments are presented to illustrate the theoretical results.