  The last two decades have seen a major expansion in the availability, size,
and precision of time-domain datasets 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 datasets. 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 modelling
choices involved in GP regression. We then review applications of GPs to
time-domain datasets 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. 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.
