This paper considers an uplink reconfigurable intelligent surface (RIS)-aided massive multiple-input multiple-output (MIMO) system, where the phase shifts of the RIS are designed relying on statistical channel state information (CSI). Considering the complex environment, the general Rician channel model is adopted for both the users-RIS links and RIS-BS links. We first derive the closed-form approximate expressions for the achievable rate which holds for arbitrary numbers of base station (BS) antennas and RIS elements. Then, we utilize the derived expressions to provide some insights, including the asymptotic rate performance, the power scaling laws, and the impacts of various system parameters on the achievable rate. We also tackle the sum-rate maximization and the minimum user rate maximization problems by optimizing the phase shifts at the RIS based on genetic algorithm (GA). Finally, extensive simulations are provided to validate the benefits by integrating RIS into conventional massive MIMO systems. Our simulations also demonstrate the feasibility of deploying large-size but low-resolution RIS in massive MIMO systems.