Emulation plays an important role in engineering design. However, most emulators such as Gaussian processes (GPs) are exclusively developed for interpolation/regression and their performance significantly deteriorates in extrapolation. To address this shortcoming, we introduce evolutionary Gaussian processes (EGPs) that aim to increase the extrapolation capabilities of GPs. An EGP differs from a GP in that its training involves automatic discovery of some free-form symbolic bases that explain the data reasonably well. In our approach, this automatic discovery is achieved via evolutionary programming (EP) which is integrated with GP modeling via maximum likelihood estimation, bootstrap sampling, and singular value decomposition. As we demonstrate via examples that include a host of analytical functions as well as an engineering problem on materials modeling, EGP can improve the performance of ordinary GPs in terms of not only extrapolation, but also interpolation/regression and numerical stability.