We present a technique for kernel-based identification of Wiener systems. We model the impulse response of the linear block with a Gaussian process. The static nonlinearity is modeled with a combination of basis functions. The coefficients of the static nonlinearity are estimated, together with the hyperparameters of the covariance function of the Gaussian process model, using an iterative algorithm based on the expectation-maximization method combined with elliptical slice sampling to sample from the posterior distribution of the impulse response given the data. The same sampling method is then used to find the posterior-mean estimate of the impulse response. We test the proposed algorithm on a benchmark of randomly-generated Wiener systems.