We present the calculation of the virtual corrections to $gg \to ZZ$ at next-to-leading order in QCD, focusing on the contribution from top quarks, which lacks a full analytic evaluation. The two-loop box diagrams are computed using a small-tansverse-momentum expansion, and the results are merged with those available in the high-energy expansion, in order to obtain an analytic description in the complete phase space.