Constructing accurate, high-dimensional molecular potential energy surfaces (PESs) for polyatomic molecules is challenging. Reproducing kernel Hilbert space (RKHS) interpolation is an efficient way to construct such PESs. However, RKHS interpolation is computationally most effective when the input energies are available on a regular grid. Thus, the number of reference energies required can become very large even for pentaatomic systems making such an approach computationally prohibitive when using high-level electronic structure calculations. Here, an efficient and robust scheme is presented to overcome these limitations and is applied to constructing high-dimensional PESs for systems with up to 10 atoms. Using energies as well as gradients reduces the number of input data required and thus keeps the number of coefficients at a manageable size. The correct implementation of permutational symmetry in the kernel products is tested and explicitly demonstrated for the highly symmetric CH4 molecule. Copyright © 2020 American Chemical Society.