This work introduces the Gaussian integration to address a smoothing problem of a nonlinear stochastic state space model. The probability densities of states at each time instant are assumed to be Gaussian, and their means and covariances are evaluated by utilizing the odd-even properties of Gaussian integral, which are further utilized to realize Rauch-Tung-Striebel (RTS) smoothing expressions. Given that the Gaussian integration provides an exact solution for the integral of a polynomial function over a Gaussian probability density function, it is anticipated to provide more accurate results than other existing Gaussian approximation-based smoothers such as extended Kalman, cubature Kalman, and unscented Kalman smoothers, especially when polynomial types of nonlinearity are present in the state space models. The developed smoothing algorithm is applied to the Van der Pol oscillator, where the nonlinearity associated with their dynamics is represented using polynomial functions. Simulation results are provided to demonstrate the superiority of the proposed algorithm.