This paper considers an uplink reconfigurable intelligent surface (RIS)-aided massive multiple-input multiple-output (MIMO) system with statistical channel state information (CSI). The RIS is deployed to help conventional massive MIMO networks serve the users in the dead zone. We consider the Rician channel model and exploit the long-time statistical CSI to design the phase shifts of the RIS, while the maximum ratio combination (MRC) technique is applied for the active beamforming at the base station (BS) relying on the instantaneous CSI. Firstly, we reveal the power scaling laws and derive the closed-form expressions for the uplink achievable rate which holds for arbitrary numbers of base station (BS) antennas. Based on the theoretical expressions, we discuss the rate performance under some special cases and provide the average asymptotic rate when using random phase shifts. Then, we consider the sum-rate maximization and the minimum user rate maximization problems by optimizing the phase shifts at the RIS. However, these two optimization problems are challenging to solve due to the complicated data rate expression. To solve these problems, we propose a novel genetic algorithm (GA) with low complexity but can achieve considerable performance. Finally, extensive simulations are provided to validate the benefits by integrating RIS into conventional massive MIMO systems. Besides, our simulations demonstrate the feasibility of deploying large-size but low-resolution RIS in massive MIMO systems.