In this paper we present a novel analysis of variance Gaussian process (ANOVA-GP) emulator for models governed by partial differential equations (PDEs) with high-dimensional random inputs. Gaussian process (GP) is a widely used surrogate modeling strategy, but it can become invalid when the inputs are high-dimensional. In this new ANOVA-GP strategy, high-dimensional inputs are decomposed into unions of local low-dimensional inputs, and principal component analysis (PCA) is applied to provide dimension reduction for each ANOVA term. We then systematically build local GP models for PCA coefficients based on ANOVA decomposition to provide an emulator for the overall high-dimensional problem. We present a general mathematical framework of ANOVA-GP, validate its accuracy and demonstrate its efficiency with numerical experiments.