The probabilistic characterization of non-Markovian responses to nonlinear dynamical systems under colored excitation is an important issue, arising in many applications. Extending the Fokker-Planck-Kolmogorov equation, governing the first-order response probability density function (pdf), to this case is a complicated task calling for special treatment. In this work, a new pdf-evolution equation is derived for the response of nonlinear dynamical systems under additive colored Gaussian noise. The derivation is based on the Stochastic Liouville equation (SLE), transformed, by means of an extended version of the Novikov-Furutsu theorem, to an exact yet non-closed equation, involving averages over the history of the functional derivatives of the non-Markovian response with respect to the excitation. The latter are calculated exactly by means of the state-transition matrix of variational, time-varying systems. Subsequently, an approximation scheme is implemented, relying on a decomposition of the state-transition matrix in its instantaneous mean value and its fluctuation around it. By a current-time approximation to the latter, we obtain our final equation, in which the effect of the instantaneous mean value of the response is maintained, rendering it nonlinear and non-local in time. Numerical results for the response pdf are provided for a bistable Duffing oscillator, under Gaussian excitation. The pdfs obtained from the solution of the novel equation and a simpler small correlation time (SCT) pdf-evolution equation are compared to Monte Carlo (MC) simulations. The novel equation outperforms the SCT equation as the excitation correlation time increases, keeping good agreement with the MC simulations.