<p>Stochastic reaction–diffusion equations play an important role in characterizing the coupled mechanisms among nonlinear reaction, spatial diffusion, and random perturbations. However, during long-time integration, the high-precision computation of the mean, variance, higher-order moments, probability density functions (PDFs), and cumulative distribution functions (CDFs) remains a challenging task. In this paper, we propose a time-splitting Fourier spectral method based on Wiener chaos expansion and Karhunen–Loève expansion for stochastic reaction–diffusion equations driven by additive noise, aiming to achieve efficient and stable uncertainty quantification over long-time integration. The computational accuracy, long-time stability, and efficiency of the proposed method are validated by two examples with analytical reference solutions, namely an OU-type stochastic benchmark and a linear stochastic heat equation. Numerical results show that the method accurately computes the mean, variance, and higher-order moments such as the third and fourth moments, and efficiently reconstructs PDFs and CDFs. Compared with the traditional Monte Carlo method, the proposed method achieves higher statistical accuracy and better probability distribution reconstruction under the same or lower computational cost over long-time integration. Furthermore, the framework is applied to one- and two-dimensional nonlinear stochastic reaction–diffusion equations to analyze the competition among noise, diffusion, and nonlinear reaction, revealing the statistical equilibrium characteristics and dynamic transition mechanisms of stochastic reaction–diffusion systems during long-time evolution.</p>