Extrinsic Gaussian process regression methods, such as wrapped Gaussian process, have been developed to analyze manifold data. However, there is a lack of intrinsic Gaussian process methods for studying complex data with manifold-valued response variables. In this paper, we first apply the parallel transport operator on Riemannian manifold to propose an intrinsic covariance structure that addresses a critical aspect of constructing a well-defined Gaussian process regression model. We then propose a novel intrinsic Gaussian process regression model for manifold-valued data, which can be applied to data situated not only on Euclidean submanifolds but also on manifolds without a natural ambient space. We establish the asymptotic properties of the proposed models, including information consistency and posterior consistency, and we also show that the posterior distribution of the regression function is invariant to the choice of orthonormal frames for the coordinate representations of the covariance function. Numerical studies, including simulation and real examples, indicate that the proposed methods work well.