In this work, we address the energy efficiency (EE) maximization problem in a downlink communication system utilizing reconfigurable intelligent surface (RIS) in a multi-user massive multiple-input multiple-output (mMIMO) setup with zero-forcing (ZF) precoding. The channel between the base station (BS) and RIS operates under a Rician fading with Rician factor K1. Since systematically optimizing the RIS phase shifts in each channel coherence time interval is challenging and burdensome, we employ the statistical channel state information (CSI)-based optimization strategy to alleviate this overhead. By treating the RIS phase shifts matrix as a constant over multiple channel coherence time intervals, we can reduce the computational complexity while maintaining an interesting performance. Based on an ergodic rate (ER) lower bound closed-form, the EE optimization problem is formulated. Such a problem is non-convex and challenging to tackle due to the coupled variables. To circumvent such an obstacle, we explore the sequential optimization approach where the power allocation vector p, the number of antennas M, and the RIS phase shifts v are separated and sequentially solved iteratively until convergence. With the help of the Lagrangian dual method, fractional programming (FP) techniques, and Lemma 1, insightful compact closed-form expressions for each of the three optimization variables are derived. Simulation results validate the effectiveness of the proposed method across different generalized channel scenarios, including non-line-of-sight (NLoS) and partially line-of-sight (LoS) conditions. This underscores its potential to significantly reduce power consumption, decrease the number of active antennas at the base station, and effectively incorporate RIS structure in mMIMO communication setup with just statistical CSI knowledge.