This paper presents a nonparametric Bayesian interpretation of kernel-based function learning with manifold regularization. We show that manifold regularization corresponds to an additional likelihood term derived from noisy observations of the function gradient along the regressors graph. The hyperparameters of the method are estimated by a suitable empirical Bayes approach. The effectiveness of the method in the context of dynamical system identification is evaluated on a simulated linear system and on an experimental switching system setup.