We calculate the two-loop QCD corrections to $gg to ZZ$ involving a closed top-quark loop. We present a new method to systematically construct linear combinations of Feynman integrals with a convergent parametric representation, where we also allow for irreducible numerators, higher powers of propagators, dimensionally shifted integrals, and subsector integrals. The amplitude is expressed in terms of such finite integrals by employing syzygies derived with linear algebra and finite field techniques. Evaluating the amplitude using numerical integration, we find agreement with previous expansions in asymptotic limits and provide ab initio results also for intermediate partonic energies and non-central scattering at higher energies.