A second-order, multilevel Monte Carlo ensemble, and hybridizable discontinuous Galerkin (MLMCE-HDG) method is proposed to solve the stochastic parabolic partial differential equations (SPDEs). By introducing an ensemble average of the diffusion coefficient, the MLMCE-HDG method results in a single discrete system with multiple right-hand-vectors, which can be solved more efficiently than a group of linear systems. A rigorous error estimate is obtained with a second-order accuracy in time and optimal convergence rate in physical space. Comparing with the multilevel Monte Carlo and hybridizable discontinuous Galerkin (MLMC-HDG) method, the MLMCE-HDG method can reduce the computational cost. Finally, we provide several numerical experiments to illustrate the theoretical results.