This paper investigates the two-timescale transmission scheme 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 nearly-passive beamforming at the RIS is adapted to the slowly-changing statistical CSI. Specifically, we first consider a system model with spatially independent Rician fading channels, which leads to tractable expressions and offers analytical insights on the power scaling laws and on the impact of various system parameters. Then, we analyze a more general system model with spatially correlated Rician fading channels and consider the impact of electromagnetic interference (EMI) caused by any uncontrollable sources present in the considered environment. For both case studies, we apply the linear minimum mean square error (LMMSE) estimator to estimate the aggregated channel from the users to the BS, utilize the low-complexity maximal ratio combining (MRC) detector, and derive a closed-form expression for a lower bound of the achievable rate. Besides, an accelerated gradient ascent-based algorithm is proposed for solving the minimum user rate maximization problem. Numerical results show that, in the considered setup, the spatially independent model without EMI is sufficiently accurate when the inter-distance of the RIS elements is sufficiently large and the EMI is mild. In the presence of spatial correlation, we show that an RIS can better tailor the wireless environment. Furthermore, it is shown that deploying an RIS in a massive MIMO network brings significant gains when the RIS is deployed close to the cell-edge users. On the other hand, the gains obtained by the users distributed over a large area are shown to be modest.