Low-thrust trajectories can benefit the search for propellant-optimal trajectories, but increases in modeling complexity and computational load remain a challenge for efficient mission design and optimization. An approach for developing models utilizing Gaussian Process (GP) regression and classification is proposed to perform computationally efficient optimization while obtaining acceptable accuracies for trajectories based on exponential sinusoid shaping. The goal of this work is to predict a combination of values of input variables which corresponds to a shape-based trajectory with the smallest total velocity increment (∆V) or propellant mass fraction (Jm). A GP classification model is constructed to assess whether a given combination of values of input variables corresponds to a feasible trajectory. GP regression models are developed to predict the total ∆V and Jm corresponding to a combination of shape parameters, which can replace the required integration along the shape. In addition, advanced regression models are developed to predict the target values while requiring only three input parameters, thereby replacing the entire shape computation. In order to develop a GP model that fits the problem at hand, the underlying functions and parameters should be selected rationally. A novel model development approach is proposed to ensure that the mean function, covariance function, likelihood function, inference method, and hyperparameters, which dominate the performance of the models, are chosen rationally in terms of mean absolute percentage error (MAPE) and prediction time. Using this approach, GP models are developed and tested on transfer trajectories from Earth to Mars and Ceres, and from Mars to Earth, and their performance, in terms of MAPE and prediction time, is compared to that of more common optimization techniques in combination with the exponential sinusoid and other shape-based methods. The results demonstrate that the computation time can significantly be reduced while achieving promising MAPE's, especially when the goal is to locate regions of feasible or near-optimal trajectories. The proposed model development procedure is tested for robustness, which provides confidence in the proposed approach. Furthermore, it is found that the models which map three input variables directly to a ∆V or Jm value perform better than the ones trained with shape information, which demonstrates the strength of GP models as applied to low-thrust trajectory optimization.