The past two decades have seen a major expansion in the availability, size, and precision of time-domain data sets in astronomy. Owing to their unique combination of flexibility, mathematical simplicity, and comparative robustness, Gaussian processes (GPs) have emerged recently as the solution of choice to model stochastic signals in such data sets. In this review, we provide a brief introduction to the emergence of GPs in astronomy, present the underlying mathematical theory, and give practical advice considering the key modeling choices involved in GP regression. We then review applications of GPs to time-domain data sets in the astrophysical literature so far, from exoplanets to active galactic nuclei, showcasing the power and flexibility of the method. We provide worked examples using simulated data, with links to the source code; discuss the problem of computational cost and scalability; and give a snapshot of the current ecosystem of open-source GP software packages. In summary: ▪ GP regression is a conceptually simple but statistically principled and powerful tool for the analysis of astronomical time series. ▪ It is already widely used in some subfields, such as exoplanets, and gaining traction in many others, such as optical transients. ▪ Driven by further algorithmic and conceptual advances, we expect that GPs will continue to be an important tool for robust and interpretable time-domain astronomy for many years to come.