In this Letter, we present for the first time a calculation of the complete next-to-leading order corrections to the $gg to ZH$ process. We use the method of small mass expansion to tackle the most challenging two-loop virtual amplitude, in which the top quark mass dependence is retained throughout the calculations. We show that our method provides reliable numeric results in all kinematic regions, and present phenomenological predictions for the total and differential cross sections at the Large Hadron Collider and its future upgrades. Our results are necessary ingredients towards reducing the theoretical uncertainties of the $pp to ZH$ cross sections down to the percent-level, and provide important theoretical inputs for future precision experimental collider programs.