The model prediction of the Gaussian process (GP) regression can be significantly biased when the data are contaminated by outliers. We propose a new robust GP regression algorithm that iteratively trims a portion of the data points with the largest deviation from the predicted mean. While the new algorithm retains the attractive properties of the standard GP as a nonparametric and flexible regression method, it can significantly reduce the influence of outliers even in some extreme cases. It is also easier to implement than previous robust GP variants that rely on approximate inference. Applied to various synthetic datasets with contaminations, the proposed method outperforms the standard GP and the popular robust GP variant with the Student's t likelihood, especially when the outlier fraction is high. Lastly, as a practical example in the astrophysical study, we show that this method can determine the main-sequence ridge line precisely in the color-magnitude diagram of star clusters.