This paper investigates the two-timescale transmission design for reconfigurable intelligent surface (RIS)-aided massive multiple-input multiple-output (MIMO) systems, where the beamforming at the base station (BS) is adapted to the rapidly-changing instantaneous channel state information (CSI), while the passive beamforming at the RIS is adapted to the slowly-changing statistical CSI. Specifically, we first propose a linear minimum mean square error (LMMSE) estimator to obtain the aggregated channel from the users to the BS in each channel coherence interval. Based on the estimated channel, we apply the low-complexity maximal ratio combining (MRC) beamforming at the BS, and then derive the ergodic achievable rate in a closed form expression. To draw design insights, we perform a detailed theoretical analysis departing from the derived ergodic achievable rate. If the BS-RIS channel is Rician distributed, we prove that the transmit power can be scaled proportionally to $1/M$, as the number of BS antennas, $M$, grows to infinity while maintaining a non-zero rate. If the BS-RIS channel is Rayleigh distributed, the transmit power can be scaled either proportionally to $1/\sqrt{M}$ as $M$ grows large, or proportionally to $1/N$ as the number of reflecting elements, $N$, grows large, while still maintaining a non-zero rate. By capitalizing on the derived expression of the data rate under the statistical knowledge of the CSI, we maximize the minimum user rate by designing the passive beamforming at the RIS. Numerical results confirm that, even in the presence of imperfect CSI, the integration of an RIS in massive MIMO systems results in promising performance gains. In addition, the obtained results reveal that it is favorable to place the RIS close to the users rather than close to the BS.