We present a novel nonparametric approach for identification of nonlinear systems. Exploiting the framework of Gaussian regression, the unknown nonlinear system is seen as a realization from a Gaussian random field. Its covariance encodes the idea of "fading" memory in the predictor and consists of a mixture of Gaussian kernels parametrized by few hyperparameters describing the interactions among past inputs and outputs. The kernel structure and the unknown hyperparameters are estimated maximizing their marginal likelihood so that the user is not required to define any part of the algorithmic architecture, e. g., the regressors and the model order. Once the kernel is estimated, the nonlinear model is obtained solving a Tikhonov-type variational problem. The Hilbert space the estimator belongs to is characterized. Benchmarks problems taken from the literature show the effectiveness of the new approach, also comparing its performance with a recently proposed algorithm based on direct weight optimization and with parametric approaches with model order estimated by AIC or BIC.